{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Survival Analysis with scikit-survival"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**scikit-survival** is a Python module for [survival analysis](https://en.wikipedia.org/wiki/Survival_analysis) built on top of [scikit-learn](https://scikit-learn.org/). It allows doing survival analysis while utilizing the power of scikit-learn, e.g., for pre-processing or doing cross-validation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "\n",
    "1. [What is Survival Analysis?](#What-is-Survival-Analysis?)\n",
    "2. [The Veterans' Administration Lung Cancer Trial](#The-Veterans'-Administration-Lung-Cancer-Trial)\n",
    "3. [Survival Data](#Survival-Data)\n",
    "4. [The Survival Function](#The-Survival-Function)\n",
    "5. [Considering other variables by stratification](#Considering-other-variables-by-stratification)\n",
    "6. [Multivariate Survival Models](#Multivariate-Survival-Models)\n",
    "7. [Measuring the Performance of Survival Models](#Measuring-the-Performance-of-Survival-Models)\n",
    "8. [Feature Selection: Which Variable is Most Predictive?](#Feature-Selection:-Which-Variable-is-Most-Predictive?)\n",
    "9. [What's next?](#What's-next?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Survival Analysis?\n",
    "\n",
    "The objective in survival analysis — also referred to as reliability analysis in engineering — is to establish a connection between covariates and the time of an event. The name *survival analysis* originates from clinical research, where predicting the time to death, i.e., survival, is often the main objective. Survival analysis is a type of regression problem (one wants to predict a continuous value), but with a twist. It differs from traditional regression by the fact that parts of the training data can only be partially observed – they are *censored*.\n",
    "\n",
    "As an example, consider a clinical study, which investigates coronary heart disease and has been carried out over a 1 year period as in the figure below.\n",
    "\n",
    "![image censoring](https://k-d-w.org/clipboard/censoring.png)\n",
    "\n",
    "Patient A was lost to follow-up after three months with no recorded cardiovascular event, patient B experienced an event four and a half months after enrollment, patient D withdrew from the study two months after enrollment, and patient E did not experience any event before the study ended. Consequently, the exact time of a cardiovascular event could only be recorded for patients B and C; their records are *uncensored*. For the remaining patients it is unknown whether they did or did not experience an event after termination of the study. The only valid information that is available for patients A, D, and E is that they were event-free up to their last follow-up. Therefore, their records are *censored*.\n",
    "\n",
    "Formally, each patient record consists of a set of covariates $x \\in \\mathbb{R}^d$ , and the time $t>0$ when an event occurred or the time $c>0$ of censoring. Since censoring and experiencing and event are mutually exclusive, it is common to define an event indicator $\\delta \\in \\{0;1\\}$ and the observable survival time $y>0$. The observable time $y$ of a right censored sample is defined as\n",
    "\n",
    "$$\n",
    "y = \\min(t, c) = \n",
    "\\begin{cases} \n",
    "t & \\text{if } \\delta = 1 , \\\\ \n",
    "c & \\text{if } \\delta = 0 .\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Consequently, survival analysis demands for models that take this unique characteristic of such a dataset into account, some of which are showcased below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Veterans' Administration Lung Cancer Trial\n",
    "\n",
    "The Veterans' Administration Lung Cancer Trial is a randomized trial of two treatment regimens for lung cancer. The [data set](http://lib.stat.cmu.edu/datasets/veteran) (Kalbfleisch J. and Prentice R, (1980) The Statistical Analysis of Failure Time Data. New York: Wiley) consists of 137 patients and 8 variables, which are described below:\n",
    "\n",
    "- `Treatment`: denotes the type of lung cancer treatment; `standard` and `test` drug.\n",
    "- `Celltype`: denotes the type of cell involved; `squamous`, `small cell`, `adeno`, `large`.\n",
    "- `Karnofsky_score`: is the Karnofsky score.\n",
    "- `Diag`: is the time since diagnosis in months.\n",
    "- `Age`: is the age in years.\n",
    "- `Prior_Therapy`: denotes any prior therapy; `none` or `yes`.\n",
    "- `Status`: denotes the status of the patient as dead or alive; `dead` or `alive`.\n",
    "- `Survival_in_days`: is the survival time in days since the treatment.\n",
    "\n",
    "Our primary interest is studying whether there are subgroups that differ in survival and whether we can predict survival times."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Survival Data\n",
    "\n",
    "As described in the section *What is Survival Analysis?* above, survival times are subject to right-censoring, therefore, we need to consider an individual's status in addition to survival time. To be fully compatible with scikit-learn, `Status` and `Survival_in_days` need to be stored as a [structured array](https://numpy.org/doc/stable/user/basics.rec.html) with the first field indicating whether the actual survival time was observed or if was censored, and the second field denoting the observed survival time, which corresponds to the time of death (if `Status == 'dead'`, $\\delta = 1$) or the last time that person was contacted (if `Status == 'alive'`, $\\delta = 0$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([( True,  72.), ( True, 411.), ( True, 228.), ( True, 126.),\n",
       "       ( True, 118.), ( True,  10.), ( True,  82.), ( True, 110.),\n",
       "       ( True, 314.), (False, 100.), ( True,  42.), ( True,   8.),\n",
       "       ( True, 144.), (False,  25.), ( True,  11.), ( True,  30.),\n",
       "       ( True, 384.), ( True,   4.), ( True,  54.), ( True,  13.),\n",
       "       (False, 123.), (False,  97.), ( True, 153.), ( True,  59.),\n",
       "       ( True, 117.), ( True,  16.), ( True, 151.), ( True,  22.),\n",
       "       ( True,  56.), ( True,  21.), ( True,  18.), ( True, 139.),\n",
       "       ( True,  20.), ( True,  31.), ( True,  52.), ( True, 287.),\n",
       "       ( True,  18.), ( True,  51.), ( True, 122.), ( True,  27.),\n",
       "       ( True,  54.), ( True,   7.), ( True,  63.), ( True, 392.),\n",
       "       ( True,  10.), ( True,   8.), ( True,  92.), ( True,  35.),\n",
       "       ( True, 117.), ( True, 132.), ( True,  12.), ( True, 162.),\n",
       "       ( True,   3.), ( True,  95.), ( True, 177.), ( True, 162.),\n",
       "       ( True, 216.), ( True, 553.), ( True, 278.), ( True,  12.),\n",
       "       ( True, 260.), ( True, 200.), ( True, 156.), (False, 182.),\n",
       "       ( True, 143.), ( True, 105.), ( True, 103.), ( True, 250.),\n",
       "       ( True, 100.), ( True, 999.), ( True, 112.), (False,  87.),\n",
       "       (False, 231.), ( True, 242.), ( True, 991.), ( True, 111.),\n",
       "       ( True,   1.), ( True, 587.), ( True, 389.), ( True,  33.),\n",
       "       ( True,  25.), ( True, 357.), ( True, 467.), ( True, 201.),\n",
       "       ( True,   1.), ( True,  30.), ( True,  44.), ( True, 283.),\n",
       "       ( True,  15.), ( True,  25.), (False, 103.), ( True,  21.),\n",
       "       ( True,  13.), ( True,  87.), ( True,   2.), ( True,  20.),\n",
       "       ( True,   7.), ( True,  24.), ( True,  99.), ( True,   8.),\n",
       "       ( True,  99.), ( True,  61.), ( True,  25.), ( True,  95.),\n",
       "       ( True,  80.), ( True,  51.), ( True,  29.), ( True,  24.),\n",
       "       ( True,  18.), (False,  83.), ( True,  31.), ( True,  51.),\n",
       "       ( True,  90.), ( True,  52.), ( True,  73.), ( True,   8.),\n",
       "       ( True,  36.), ( True,  48.), ( True,   7.), ( True, 140.),\n",
       "       ( True, 186.), ( True,  84.), ( True,  19.), ( True,  45.),\n",
       "       ( True,  80.), ( True,  52.), ( True, 164.), ( True,  19.),\n",
       "       ( True,  53.), ( True,  15.), ( True,  43.), ( True, 340.),\n",
       "       ( True, 133.), ( True, 111.), ( True, 231.), ( True, 378.),\n",
       "       ( True,  49.)],\n",
       "      dtype=[('Status', '?'), ('Survival_in_days', '<f8')])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.datasets import load_veterans_lung_cancer\n",
    "\n",
    "data_x, data_y = load_veterans_lung_cancer()\n",
    "data_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily see that only a few survival times are right-censored (`Status` is `False`), i.e., most veteran's died during the study period (`Status` is `True`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Survival Function\n",
    "\n",
    "A key quantity in survival analysis is the so-called survival function, which relates time to the probability of surviving beyond a given time point.\n",
    "\n",
    "> Let $T$ denote a continuous non-negative random variable corresponding to a patient’s survival time. The survival function $S(t)$ returns the probability of survival beyond time $t$ and is defined as\n",
    "> $$ S(t) = P (T > t). $$\n",
    "\n",
    "If we observed the exact survival time of all subjects, i.e., everyone died before the study ended, the survival function at time $t$ can simply be estimated by the ratio of patients surviving beyond time $t$ and the total number of patients:\n",
    "\n",
    "$$\n",
    "\\hat{S}(t) = \\frac{ \\text{number of patients surviving beyond $t$} }{ \\text{total number of patients} }\n",
    "$$\n",
    "\n",
    "In the presence of censoring, this estimator cannot be used, because the numerator is not always defined. For instance, consider the following set of patients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Status</th>\n",
       "      <th>Survival_in_days</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>True</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>True</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>True</td>\n",
       "      <td>20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>False</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>True</td>\n",
       "      <td>59.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Status  Survival_in_days\n",
       "1    True               8.0\n",
       "2    True              10.0\n",
       "3    True              20.0\n",
       "4   False              25.0\n",
       "5    True              59.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "pd.DataFrame.from_records(data_y[[11, 5, 32, 13, 23]], index=range(1, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the formula from above, we can compute $\\hat{S}(t=11) = \\frac{3}{5}$, but not $\\hat{S}(t=30)$, because we don't know whether the 4th patient is still alive at $t = 30$, all we know is that when we last checked at $t = 25$, the patient was still alive.\n",
    "\n",
    "An estimator, similar to the one above, that *is* valid if survival times are right-censored is the [Kaplan-Meier estimator](https://en.wikipedia.org/wiki/Kaplan%E2%80%93Meier_estimator)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'time $t$')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlsUlEQVR4nO3de1SUdf4H8PcImsniBVQExhuO4oCgIqSImcoKRi4qti6meeUQalnrVtp2ftaWm7qdXXPVZMmWNFPyuJvsmmCG1pFQLspqB5HlIBqDFhdRSQMRvr8/bGa5DjwPzzC39+sczuGZeXjm8+105uP39vmqhBACREREbehm7gCIiMiyMVEQEZFRTBRERGQUEwURERnFREFEREYxURARkVFmTRQrVqzAwIEDMWbMmFbfF0Jg7dq10Gg08Pf3x/nz57s4QiIiMmuiWLZsGVJTU9t8PyUlBYWFhSgsLERCQgJWrVrVhdERERFg5kQxdepUuLi4tPl+cnIylixZApVKhUmTJuHWrVu4ceNGF0ZIREQWPUdRWlqKwYMHG67VajVKS0vNGBERkf1xNHcAxrRWXUSlUrV4LSEhAQkJCQCAy5cvY/To0SaPjYjIlly9ehUVFRWtvmfRiUKtVqOkpMRwrdPp4OHh0eK+2NhYxMbGAgACAwORk5PTZTESEdmCwMDANt+z6KGnyMhI7Nu3D0IInD17Fn369IG7u7vJPu8P/85D3P5zJns+EZE1MmuPYuHChfjqq69QUVEBtVqNP/zhD6irqwMAxMXFISIiAseOHYNGo0GvXr2QmJho0nguXb+DOzV1Jv0MIiJrY9ZEcfDgQaPvq1Qq7Nq1q4uiISKi1lj00BMREZkfE0Uz+TeqcSDzO3OHQURkMZgoGpkzzhMA8PHZq+YNhIjIgjBRNPLMxCHQujubOwwiIovCREFEREYxURARkVFMFEREZBQTRSu48omI6H+YKJoJGdEfAJD8H1apJSICmChaCNW6QevuzFIeREQ/Y6IgIiKjmCiIiMgoJopmgke4mjsEIiKLwkRBRERGMVEQEZFRTBRGnCmqNHcIRERmx0RBRERGMVG0If9GNdLyfzB3GEREZsdE0Qr97uxviirMHAkRkfkxUbRCvzsb4DwFERETRSu4l4KI6H+YKIiIyCgmCiIiMoqJgoiIjGKiICIio5gojOBeCiIiJoo2cS8FEdFDTBRtaLyXgojInjFREBGRUUwUHcDd2URkz5go2sEJbSKyd51KFHfv3kV9fb1SsVgcTmgTEQGOUm5uaGhAUlISPvnkE2RnZ+ORRx5BbW0tBgwYgIiICMTGxmLkyJGmirVL6es9MUkQkb2T1KOYPn06ioqKsHnzZnz//fcoKSlBWVkZTp8+jUmTJmHDhg3Yv3+/qWI1K85TEJG9ktSj+PLLL9G9e3dcu3YN3br9L8e4uLhg/vz5mD9/Purq6hQPkoiIzEdSj6J79+4AgHnz5rV47+zZs03uISIi2yApURw6dAgbNmxAdXU18vPzm0xkx8bGSv7w1NRUeHt7Q6PRYMuWLS3ev337Nn71q19h7Nix8PX1RWJiouTPICKizpE09BQSEoKamhrs2bMH69atQ0FBAfr27QsPDw88+uijkj64vr4ea9aswYkTJ6BWqxEUFITIyEj4+PgY7tm1axd8fHzw73//G+Xl5fD29saiRYvQo0cPSZ9FRETySUoUnp6eWLJkCUaMGIGQkBAAwM2bN1FcXIzRo0dL+uCsrCxoNBp4eXkBAKKjo5GcnNwkUahUKlRXV0MIgR9//BEuLi5wdJQUMhERdZKkb10hBFQqlSFJAA8nsl1cXFrc057S0lIMHjzYcK1Wq5GZmdnknueffx6RkZHw8PBAdXU1Pv300yaT6HoJCQlISEgAAJSXl0tpEhERtUPy8tgdO3bgu+++a/L6/fv3cfLkSSxduhR79+7t0LOEEC1ea55gjh8/jnHjxuH69ev4z3/+g+effx537txp8XexsbHIyclBTk4OBgwYIKFF0nCJLBHZI0mJIjU1FQ4ODli4cCE8PDzg4+OD4cOHY+TIkTh48CB++9vfYtmyZR16llqtRklJieFap9PBw8OjyT2JiYmIioqCSqWCRqPB8OHDcfnyZSkhd4p+0x3LeBCRPZM09NSzZ0+sXr0aq1evRl1dHSoqKvDoo4+ib9++kj84KCgIhYWFKC4uhqenJ5KSknDgwIEm9wwZMgRpaWl4/PHH8cMPP6CgoMAwp9FVQkb0R/6NanxTVIFQrVuXfjYRkSUw2qO4dOkSFi9ebLgODQ1FXl4egIf7JbKzs7Fz505kZWVJ/mBHR0fs3LkT4eHh0Gq1WLBgAXx9fREfH4/4+HgAwP/93/8hIyMDfn5+CA0NxdatW9G/f3/Jn9UZzc+l4PATEdkboz2K0NBQnDlzxnCt0+ng6+sLAMjIyMCzzz6L3/zmN1i2bBn++Mc/troRz5iIiAhEREQ0eS0uLs7wu4eHB7744gtJzyQiImUZ7VF88cUXeP311w3XvXv3Nvy+b98+xMXFISEhAV999RW2bt1quiiJiMhsjCYKPz8/fPLJJ4ZrjUaDw4cPo6ysDEeOHMGcOXMAAAMHDkRtba1pIzUzTmgTkb2StOpp27Zt+Nvf/gZPT08EBARg8uTJAIC6ujr8+OOPJgnQ3IJHuPJcCiKya5JWPQ0aNAgnTpxAQ0NDk41vp06dwvTp0xUPzlKEat2YJIjIbsmqh9F8d3RYWBjCwsIUCYiIiCwLz8wmIiKjmChk4F4KIrInTBRERGSUpDkKZ2fnVivD6ivGtlawj4iIrJukRFFdXW2qOKyCfi9FqNbNMPykLxxIRGSrZJ8CVFVVhcLCQtTU1Bhemzp1qiJBWSIWByQieyUrUezZswfbt2+HTqfDuHHjcPbsWQQHB+PkyZNKx2cxuJeCiOyVrMns7du3Izs7G0OHDsWpU6eQm5tr0gODLEnzUh5cAUVEtk5WoujZsyd69uwJAKitrcXo0aNRUFCgaGCWiKU8iMgeyRp6UqvVuHXrFubOnYuZM2eiX79+LU6nsyXBI1xxpqiSw09EZJdkJYrPPvsMAPDmm29i+vTpuH37NmbNmqVoYNbkTFElVz8Rkc2SNfS0bds26HQ6AMATTzyByMhI9OjRQ9HAiIjIMshKFHfu3EF4eDgef/xx7Nq1Cz/8wHMaiIhslaxE8cYbbyAvLw+7du3C9evX8cQTT+CXv/yl0rFZFa5+IiJb1alaTwMHDsSgQYPg6uqKsrIypWIiIiILIitR7N69G9OmTUNoaCgqKirwwQcf4OLFi0rHZlEaT1bzWFQisieyVj1du3YN7733HsaNG6dwOJaPpTyIyN7IShRbtmxROg6rwb0URGRvJA09TZkyBcDDcuO9e/c2/OiviYjI9kjqUaSnpwNguXEAuFZ5D28dzUPIiP4cgiIimyZ7w11paanSsViNkBH9MdS1F65V3uMwFBHZPNkb7sLCwux2w12o1g0bZ/tiqGsvc4dCRGRy3HDXSY2XynLTHRHZIm64k6B54T+WHScie8ANd50QqnWD1t3Z3GEQEZmU5H0UQgjk5OTY7Ya71uiHn7j6iYhskeQehUqlQm5uLpPEzzj8RES2TtbQU3BwMLKzs5WOxSo0n6fg8BMR2TpZJTxOnTqF+Ph4DBs2DE5OThBCQKVS2eU8hV7j1U887Y6IbImsRJGSkqJ0HFaNhQKJyJbJShR79+5t9fWNGzdKek5qaipefPFF1NfXIyYmBhs2bGhxz1dffYWXXnoJdXV16N+/P77++ms5IZtU80KBPEObiGyJrETh5ORk+L2mpgZHjx6FVquV9Iz6+nqsWbMGJ06cgFqtRlBQECIjI+Hj42O459atW1i9ejVSU1MxZMgQq9qrwWRBRLZCVqL43e9+1+T65ZdfRmRkpKRnZGVlQaPRwMvLCwAQHR2N5OTkJoniwIEDiIqKwpAhQwA83OBHRERdq1M7s/Xu3buHK1euSPqb0tJSDB482HCtVqtbFBr873//i6qqKkybNg0TJkzAvn37lAi3y5wpqmRZDyKyerJ6FH5+flCpVAAeDiGVl5dLnp8QQrR4Tf9MvQcPHuDcuXNIS0vDTz/9hODgYEyaNAmjRo1qcl9CQgISEhIAAOXl5ZLiICIi42QliqNHj/7vAY6OcHNzg6OjtEep1WqUlJQYrnU6HTw8PFrc079/fzg5OcHJyQlTp07FhQsXWiSK2NhYxMbGAgACAwOlNkey4BGuknoKnK8gImsma+gpKysLLi4uGDp0KBITE7FgwQKcP39e0jOCgoJQWFiI4uJi3L9/H0lJSS3mOebMmYPTp0/jwYMHuHfvHjIzMyVPmlsKDkMRkbWSlSjefvttODs7Iz09HcePH8fSpUuxatUqSc9wdHTEzp07ER4eDq1WiwULFsDX1xfx8fGIj48HAGi1WsyaNQv+/v547LHHEBMTgzFjxsgJmYiIZFKJ1iYL2jF+/Hjk5ubitddeg5+fH5555hnDa+YWGBiInJwck39O897BW0fzkH+jGjFThhvddMchKCKyRMa+O2X1KDw9PfHcc8/h0KFDiIiIQG1tLRoaGjoVpLUJHuHa5EufxQGJyFbJShSHDh1CeHg4UlNT0bdvX9y8eRPvvvuu0rFZFRYHJCJbJWvVU69evRAVFWW4dnd3h7u7u2JB2TKugCIia6PIhjt7xi99IrJ1TBRmwGWyRGRNJCWKZ599FgCwfft2kwRjT5gsiMhaSEoU586dw7Vr1/D3v/8dVVVVuHnzZpMfe9V4+KnxAUZERLZA0mR2XFwcZs2ahStXrmDChAlN6jWpVCrJhQFtjf4Aoz3pxQDAQ4yIyCZI6lGsXbsW+fn5WLFiBa5cuYLi4mLDj70nCeBhYoiZMhwA91MQke2QtTx29+7duHDhAk6fPg0AmDp1Kvz9/RUNzFrpT7vTD0EZ61VwqSwRWQNZq57++te/YtGiRSgrK0NZWRkWLVqEHTt2KB2bVeEubSKyVbJqPfn7++PMmTOGI1Hv3r2L4OBgXLx4UfEApeqqWk+tabyS6a2jebhWeQ9DXXshZET/NnsW7FEQkSUw9t0pa+hJCAEHBwfDtYODQ6sHEdmzh72Kh0NQ+TeqAbQ+ua1PLkwYRGSpZCWK5cuXY+LEiZg3bx4A4MiRI1i5cqWigVm7UK0bQrVuSMv/AXvSi/FNUQVXQRGRVZKVKNatW4dp06YhPT0dQggkJiZi/PjxSsdmE/ST29cq7+Gto3lGh6GIiCyRrEQBAAEBAQgICFAyFpulH4a6VnkPAHsWRGRdWOupC4Rq3bBxti+GuvYydyhERJIxUVgI1n4iIkslK1Hs3LkTVVVVSsdCREQWSFai+P777xEUFIQFCxYgNTWVS2MlYNFAIrI2shLFpk2bUFhYiJUrV+Kjjz7CyJEj8fvf/x5FRUVKx2dV2tsLwR3bRGSNZM9RqFQqDBo0CIMGDYKjoyOqqqrw9NNP49VXX1UyPpuiP1e7rV7FmaJKzlUQkcWRXetpwoQJePXVVxESEoJvv/0Wu3fvxrlz5/CPf/xD6Rhtir5XsSe9mENQRGQVZO2jqKiowD//+U8MHTq0yevdunXD0aNHFQnMVun3UHC3NhFZC1k9itra2hZJYv369QAArVbb+ahsnH4IiojIGshKFCdOnGjxWkpKSqeDISIiyyMpUezevRt+fn4oKCiAv7+/4Wf48OE8uOhnUqrAGpvUJiKyFJLmKJ555hk8+eSTeO2117BlyxbD687OznBxcVE8OFumP1+b8xREZOkkJYo+ffqgT58+OHjwoKnisQnBI1zb7RXoq8q2hcekEpGlkDT0NGXKFAAPexC9e/c2/Oiv6X86+iWvLz/OpbJEZKkk9SjS09MBANXV1SYJxt6w/DgRWQNWjzWh4BGuRnsWLD9ORNZAUo/C2dkZKpWq1SKAKpUKd+7cUSwwIiKyDJISBYec5OnI5HZr9H/DSW0iMidFJrP1P9S29oahWH6ciCyVpETReDL7zp07LX6kSk1Nhbe3NzQaTZN9Gc1lZ2fDwcEBhw8flvwZ1qC98uPcgEdE5mS2yez6+nqsWbMGKSkpuHTpEg4ePIhLly61et/69esRHh5uhii7Bms/EZElk5Uoampq8Je//AVRUVGYP38+tm3bhpqaGknPyMrKgkajgZeXF3r06IHo6GgkJye3uG/Hjh2YP38+Bg4cKCdUm8GzKojIXGQliiVLliAvLw8vvPACnn/+eeTn5+PZZ5+V9IzS0lIMHjzYcK1Wq1FaWtrins8++wxxcXFywrQ63HxHRJZI1nkUBQUFuHDhguF6+vTpGDt2rKRntLXEtrGXXnoJW7duhYODg9FnJSQkICEhAQBQXl4uKY6u1tYKKG6+IyJLJStRjB8/HmfPnsWkSZMAAJmZmQgJCZH0DLVajZKSEsO1TqeDh4dHk3tycnIQHR0N4OFhSceOHYOjoyPmzp3b5L7Y2FjExsYCAAIDA6U2xyKEat0QqnXDW0fzzB0KEVETkhKFn58fVCoV6urqsG/fPgwZMgQA8N1338HHx0fSBwcFBaGwsBDFxcXw9PREUlISDhw40OSe4uJiw+/Lli3D7NmzWyQJayR3XwXAvRVE1PUkJQoljzl1dHTEzp07ER4ejvr6eqxYsQK+vr6Ij48HAJufl+hMsiAi6kqSEkXj40+rqqpQWFjYZLVT8+NR2xMREYGIiIgmr7WVID766CNJz7Zm+s13nKcgIksga9XTnj17MHXqVISHh+ONN95AeHg43nzzTYVDs0/tbb7TY2+EiLqKrESxfft2ZGdnY+jQoTh16hRyc3MxYMAApWOzea3NM+g337GkBxFZClmJomfPnujZsycAoLa2FqNHj0ZBQYGigdmzjvYqiIi6gqzlsWq1Grdu3cLcuXMxc+ZM9OvXr8XSVpKvvWNSiYi6kqxE8dlnnwEA3nzzTUyfPh23b9/GrFmzFA2M2selskTUFWQlipqaGrz//vtIT0+HSqXClClT0NDQoHRsdkH/Jc/JaSKyVGar9UTt62jtJyYZIjIls9V6IuNY+4mILIWsHoW+1pOenFpPZFyo1g0bZ/tiqGuvDi2VZa+CiEzFbLWeqGNCRvRH/o1qfFPEXgURmYfZaj1RU23VftIvldXPV4SM6N9mwuAqKCIyBdm1ni5cuIDTp08DAB5//HHOUZgQ5yuIyJxkl/BYtGgRysrKUFZWhsWLF2PHjh1Kx0Y/azxfQUTU1WStevrwww+RmZkJJycnAMD69esRHByMF154QdHg7E1H9lR0ZAiKiEhJshKFEKLJ8aQODg6tHm1KyuIQFBGZg6xEsXz5ckycOBHz5s0DABw5cgQrV65UNDBqicelEpE5yEoU69atw7Rp05Ceng4hBBITEzF+/HilYyOZmg9dcRUUEXWG5EQhhIBOp0NAQAACAgJMERN1AE/BI6KuInnVk0qlwty5c00QCum11wOQel7FmaJKww8RkVSylsdOmjQJ2dnZSsdCHaQ/BY+IqCvImqM4deoU4uPjMWzYMDg5OUEIAZVKhYsXLyodn91qa6d2Z3H+goikkpUoUlJSlI6DZOCeCiLqCrIShZubW4uDi1atWqV0bGSEUnsqzhRVsldBREbx4CIrJbUMORGRXDy4yMrpy5DvSS8GAFk9C85bEJExPLjIggWPcG33SztU64aYKcMBdHy5LBGRFCoho0iTVqtFQUFBk4OLtFotunXrZvbVT4GBgcjJyTHb55tKeyug3jqah2uV9zDUtZfik9vsYRDZPmPfnbKGnlJTUzsVECmPBQOJyFRkJYrGBxhR12ivBLkpCwZyZRSRfZM1R0FERPZDVo+C7E/jngx7F0T2RbEexffff6/Uo4iIyIIo1qNYuXIlPv/8c6UeR22QclyqntKroNi7ILIviiUKJgnLoF/9pJd/oxr5N6oNeyxYF4qIpJI19LR+/foOvUZdT1/aQ/8TM2W4oST5tcp7im/K41kXRLZP1oa7gIAAnD9/vslr/v7+FlFm3FY33LVG6pdz4015ekr2MDgMRWS9jH13SupR7N69G35+frh8+TL8/f0NP8OHD4efn5/kwFJTU+Ht7Q2NRoMtW7a0eP+TTz4xfMbkyZOb1Jci6UJG9G+SJEzRwyAi2yOpR3H79m1UVVXhtddea/LF7uzsDBcXF0kfXF9fj1GjRuHEiRNQq9UICgrCwYMH4ePjY7gnIyMDWq0W/fr1Q0pKCt58801kZmYafS57FB331tE85N+oRsyU4Sabt2Avg8g6KNaj6NOnD4YNG4aoqCi4uLhg6NCh+PjjjxETE4Pc3FxJQWVlZUGj0cDLyws9evRAdHQ0kpOTm9wzefJk9OvXD8DD41d1Op2kz7B1HSkaaIzUs7eJyD7Jmsx+++234ezsjPT0dBw/fhxLly5FXFycpGeUlpZi8ODBhmu1Wo3S0tI27//www/x5JNPtvpeQkICAgMDERgYiPLycklx2DP92dv65bT6HyXPtuBEN5H1k5UoHBwcADxcErtq1SrMmTMH9+/fl/SM1ka8VCpVq/eeOnUKH374IbZu3drq+7GxscjJyUFOTg4GDBggKQ5b0NleBectiMgYWfsoPD098dxzz+HEiRNYv349amtr0dDQIOkZarUaJSUlhmudTgcPD48W9128eBExMTFISUmBqyvHu9sSPMJV1r/c9cUE9fQro5TesMfCgkTWS1aP4tChQwgPD8fx48fRt29f3Lx5E++++66kZwQFBaGwsBDFxcW4f/8+kpKSEBkZ2eSe7777DlFRUfj4448xatQoOaGSROxhEFFzsnoUjz76KO7evYuDBw9i48aNqKurQ9++faV9sKMjdu7cifDwcNTX12PFihXw9fVFfHw8ACAuLg5vvfUWKisrsXr1asPf2MuKJjk6Ut6jPa31MPRncivRqzCGPQ4iyyRrw92qVavQrVs3nDx5Evn5+aiqqkJYWBiys7NNEaMk9rQ8ti1KTh6n5f+APenF0Lo7Y+NsX8We2xomCiLzUfyEu8zMTJw/fx7jx48HAPTr10/yZDZZh1CtG74pqmgyb8F6UUT2RVai6N69O+rr6w2rlMrLy9GtG89AslWNCw2a8qhVDk0RWSZZiWLt2rWYN28eysrK8Prrr+Pw4cPYtGmT0rGRhWg8b9F8VRR7F0S2T1aiWLRoESZMmIC0tDQIIXDkyBFotVqlYyOZmv/LW8k5i67qXbSmtXawl0FkerImsy0ZJ7ONUzJpNK9Ga47eBRMFkTIUn8wmAszbu9AzlviYRIiUwURhZ5QcljI2d6HHOQwi68dEQYpofgQrYL5eBhEpi3MUZLLqro3nMMzVs+DwE1HHcI6CjDLVKil9L4M9CyLrxl1yZDKhWjdsnO2Loa69DPWiuhrPwyDqPPYoqAWlexghI/oj/0Y1vilir4LIGjFRkMm1Vi+qNVwhRWSZOPRE7VJiQrj5ORfNmfrcCw4/EcnHHgV1SONkocRJes21tQ+jMfY4iMyDiYIsQmv7MBpTYuVUWwmOS2iJjGOiIMk627tojdweB3sZRKbHREFWgTu/icyHO7NJMV09Ydy8em1zXdnb4PAVWTvuzKYu0dqXpSmTh7F5DfY2iJTDREFWy9i8xltH8wy7wZksiDqHiYJMSt/L6Ophqa7eDd5V7eMQF5kDN9xRl+jqL7hQrRu07s6GlVLmqDNFZCvYoyCbZYvVazvac2HPg5TEREFdpr0vL6WHb/RzGJyvIOocJgqyGKZKJKxeS9Q5TBRk89qrXmuLu7tZBNG+mHqokYmCrEZn9mm0tefCluYviEyFiYKsmrF/STVOIm3tuehI1dqOsMVeCZEeEwXZtfaq1nYEeyVkbvp/FJlqCIqJgmxW8AjXdoem2qta2xFK9Uo6gj0XMgcmCrJpXbEzXIleSUew50LmwkRB1ElK9Eo6QkrPhT0PUhJLeJBdsIWdyu2dO65n6vPHyf6wR0F2o61kYS17Djrac+nKOROyDKbuQZo1UaSmpuLFF19EfX09YmJisGHDhibvCyHw4osv4tixY+jVqxc++ugjBAQEmClaslVK9DYsKdl01ZwJWYaumLsyW6Kor6/HmjVrcOLECajVagQFBSEyMhI+Pj6Ge1JSUlBYWIjCwkJkZmZi1apVyMzMNFfIRG3qqqGtjiSkrpozIcvQuAf57KRheGbiEMU/w2yJIisrCxqNBl5eXgCA6OhoJCcnN0kUycnJWLJkCVQqFSZNmoRbt27hxo0bcHd3N1fYRGbV0YRkST0cMq3GVZKT/1NqW4mitLQUgwcPNlyr1eoWvYXW7iktLWWiIGqHLUzeU8fo9wuZck7KbIlCCNHiNZVKJfkeAEhISEBCQgIA4PLlywgMDJQVU3l5OQYMGCDrb60V22wf2Gb7cKm8HIEfyGvz1atX23zPbIlCrVajpKTEcK3T6eDh4SH5HgCIjY1FbGxsp2MKDAxETk5Op59jTdhm+8A22wdTtdls+yiCgoJQWFiI4uJi3L9/H0lJSYiMjGxyT2RkJPbt2wchBM6ePYs+ffpw2ImIqIuZrUfh6OiInTt3Ijw8HPX19VixYgV8fX0RHx8PAIiLi0NERASOHTsGjUaDXr16ITEx0VzhEhHZLbPuo4iIiEBEREST1+Li4gy/q1Qq7Nq1q8viUWL4ytqwzfaBbbYPpmqzSrQ2Y0xERPQz1noiIiKjmCh+lpqaCm9vb2g0GmzZssXc4SimpKQE06dPh1arha+vL7Zv3w4AuHnzJmbOnImRI0di5syZqKqqMvzN5s2bodFo4O3tjePHj5sr9E6pr6/H+PHjMXv2bAC2314AuHXrFp5++mmMHj0aWq0WZ86csel2b9u2Db6+vhgzZgwWLlyImpoam2zvihUrMHDgQIwZM8bwmpx2njt3Dn5+ftBoNFi7dm2r2w/aJEg8ePBAeHl5iaKiIlFbWyv8/f1FXl6eucNSxPXr18W5c+eEEELcuXNHjBw5UuTl5YlXXnlFbN68WQghxObNm8Wrr74qhBAiLy9P+Pv7i5qaGnHlyhXh5eUlHjx4YLb45frzn/8sFi5cKJ566ikhhLD59gohxJIlS8QHH3wghBCitrZWVFVV2Wy7dTqdGDZsmLh3754QQohf//rXIjEx0Sbb+/XXX4tz584JX19fw2ty2hkUFCQyMjJEQ0ODmDVrljh27FiHY2CiEEJkZGSIsLAww/U777wj3nnnHTNGZDqRkZHiiy++EKNGjRLXr18XQjxMJqNGjRJCtGx7WFiYyMjIMEuscpWUlIgZM2aItLQ0Q6Kw5fYKIcTt27fFsGHDRENDQ5PXbbXdOp1OqNVqUVlZKerq6sRTTz0ljh8/brPtLS4ubpIopLbz+vXrwtvb2/D6gQMHRGxsbIc/n0NPaLtUiK25evUqcnNzMXHiRPzwww+GPSnu7u4oKysDYBv/LV566SX86U9/Qrdu//vf25bbCwBXrlzBgAEDsHz5cowfPx4xMTG4e/euzbbb09MTL7/8MoYMGQJ3d3f06dMHYWFhNtve5qS2s7S0FGq1usXrHcVEgY6XCrFmP/74I+bPn4/33nsPvXv3bvM+a/9vcfToUQwcOBATJkzo0P3W3l69Bw8e4Pz581i1ahVyc3Ph5ORkdK7N2ttdVVWF5ORkFBcX4/r167h79y7279/f5v3W3t6OaqudnW0/EwU6XirEWtXV1WH+/PlYtGgRoqKiAABubm64ceMGAODGjRsYOHAgAOv/b/HNN9/gX//6F4YNG4bo6GicPHkSixcvttn26qnVaqjVakycOBEA8PTTT+P8+fM22+4vv/wSw4cPx4ABA9C9e3dERUUhIyPDZtvbnNR2qtVq6HS6Fq93FBMFOlZOxFoJIbBy5UpotVqsW7fO8HpkZCT27t0LANi7dy/mzJljeD0pKQm1tbUoLi5GYWEhHnvsMbPELsfmzZuh0+lw9epVJCUlYcaMGdi/f7/Ntldv0KBBGDx4MAoKCgAAaWlp8PHxsdl2DxkyBGfPnsW9e/cghEBaWhq0Wq3Ntrc5qe10d3eHs7Mzzp49CyEE9u3bZ/ibDpE3tWJ7Pv/8czFy5Ejh5eUlNm3aZO5wFHP69GkBQPj5+YmxY8eKsWPHis8//1xUVFSIGTNmCI1GI2bMmCEqKysNf7Np0ybh5eUlRo0aJWllhKU5deqUYTLbHtqbm5srJkyYIPz8/MScOXPEzZs3bbrdGzduFN7e3sLX11csXrxY1NTU2GR7o6OjxaBBg4Sjo6Pw9PQUe/bskdXO7Oxs4evrK7y8vMSaNWtaLHwwhjuziYjIKA49ERGRUUwURERkFBMFEREZxURBRERGMVEQEZFRTBRERGQUEwURERnFREEk0a1bt/D+++8bridPnmySz9HpdPj0009N8mwiKZgoiCRqnigyMjJM8jlpaWk4f/68SZ5NJAUTBZFEGzZsQFFREcaNG4dXXnkFv/jFLwA8LOM+evRoxMTEYMyYMVi0aBG+/PJLhISEYOTIkcjKyjI8Y//+/Xjssccwbtw4PPfcc6ivr2/yGenp6Vi3bh0OHz6McePGobi4uEvbSNSEgiVJiOxC80NknJycDK87ODiIixcvivr6ehEQECCWL18uGhoaxJEjR8ScOXOEEEJcunRJzJ49W9y/f18IIcSqVavE3r17W3xOeHi4+Pbbb03fIKJ2OJo7URHZkuHDh8PPzw8A4Ovri9DQUKhUKvj5+eHq1asAHg4pnTt3DkFBQQCAn376yVAmurGCggJ4e3t3WexEbWGiIFLQI488Yvi9W7duhutu3brhwYMHAB6Wfl+6dCk2b97c5nMqKyvRp08fdO/e3bQBE3UA5yiIJHJ2dkZ1dbXsvw8NDcXhw4cNx1fevHkT165da3JPcXGxVR+sQ7aFiYJIIldXV4SEhGDMmDF45ZVXJP+9j48PNm3ahLCwMPj7+2PmzJmG08r0Ro8ejYqKCowZM8Zkq6qIOornURARkVHsURARkVFMFEREZBQTBRERGcVEQURERjFREBGRUUwURERkFBMFEREZxURBRERG/T8EDYr0dI50RAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sksurv.nonparametric import kaplan_meier_estimator\n",
    "\n",
    "time, survival_prob, conf_int = kaplan_meier_estimator(\n",
    "    data_y[\"Status\"], data_y[\"Survival_in_days\"], conf_type=\"log-log\"\n",
    ")\n",
    "plt.step(time, survival_prob, where=\"post\")\n",
    "plt.fill_between(time, conf_int[0], conf_int[1], alpha=0.25, step=\"post\")\n",
    "plt.ylim(0, 1)\n",
    "plt.ylabel(r\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated curve is a step function, with steps occurring at time points where one or more patients died. From the plot we can see that most patients died in the first 200 days, as indicated by the steep slope of the estimated survival function in the first 200 days."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Considering other variables by stratification\n",
    "\n",
    "### Survival functions by treatment\n",
    "\n",
    "Patients enrolled in the Veterans' Administration Lung Cancer Trial were randomized to one of two treatments: `standard` and a new `test` drug. Next, let's have a look at how many patients underwent the standard treatment and how many received the new drug."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Treatment\n",
       "standard    69\n",
       "test        68\n",
       "Name: count, dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_x[\"Treatment\"].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roughly half the patients received the alternative treatment.\n",
    "\n",
    "The obvious questions to ask is:\n",
    "> *Is there any difference in survival between the two treatment groups?*\n",
    "\n",
    "As a first attempt, we can estimate the survival function in both treatment groups separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x120beee40>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz+klEQVR4nO3deVRUV7o28KcEDEojihNCKYOoDBaigorGKCHOinFKTDTOl4hGb+JNWnN7pTN5o4nd1zbRDk1iG+1Evcavo9EEjEG0pREVNZoGNLQCsXAGRByY9/cHqQpDUXAOp6jp+a3Faqo4dc570i5ezn73frdKCCFARETUiDbmDoCIiCwbEwURERnFREFEREYxURARkVFMFEREZBQTBRERGWXWRLFo0SJ069YN/fv3N/hzIQRWrlwJf39/hISE4OzZs60cIRERmTVRLFiwAImJiY3+PCEhAdnZ2cjOzkZ8fDxiY2NbMToiIgLMnCieeOIJuLu7N/rz/fv3Y968eVCpVBg2bBju3r2L69evt2KERERk0TWK/Px89OzZU/9arVYjPz/fjBEREdkfR3MHYIyh7iIqlarBe/Hx8YiPjwcAXLx4EQEBAbKv+aCsEi6PWfR/FiIixeXm5uLOnTsGf2bRvxHVajWuXr2qf63VauHp6dnguJiYGMTExAAAwsLCkJ6eLvuaJy4XIKJ3Z9mfJyKyRmFhYY3+zKKHnqKjo7Fjxw4IIZCWlgY3Nzf06NHDZNd7+0AG3jmYgZ0nfzbZNYiIrI1Znyiee+45HD16FHfu3IFarcbbb7+NiooKAMDSpUsxceJEfPvtt/D390f79u2xbds2k8eUdb0Ef0vLxfNDe5n8WkRE1sCsiWLXrl1Gf65SqbBly5ZWigZ4c0ow0i/morLVrkhEZPksukZBRA1VVFRAq9WitLTU3KGQFXJ2doZarYaTk1OzP8NEUVvCGiwpy0Fc2+XmjoSoUVqtFq6urvDx8TE4C5CoMUIIFBQUQKvVwtfXt9mfs+hidqu78SP8qvPgWF5i7kiIGlVaWorOnTszSZBkKpUKnTt3lvw0ykRBZIWYJEguOf92mCjq0VRnof/dJE6RJWpEQUEBQkNDERoaCg8PD3h5eelfl5eXyzrn0aNHkZqaqnCkrX+dzz77DNeuXVPsfKNHj27RurCjR49i8uTJLY6DiaI2zUwAwFSHVOz/ga1CiAzp3LkzfvjhB/zwww9YunQpXnnlFf3rtm3borJS+rxBJgplVFVVmeS8LGbXFrYQD/75F7gUAw/vFZo7GiKrsWDBAri7u+PcuXMYNGgQli1bhuXLl+P27dto3749PvnkEwQEBODAgQNYu3YtysvL0blzZ3zxxRd49OgR4uLi4ODggM8//xwfffQRtm7dinbt2uHixYvIy8vDtm3bsH37dpw4cQJDhw7FZ599BgD47rvv8Oabb6KsrAy9e/fGtm3b8Jvf/AY+Pj6YP38+Dhw4gIqKCnz55ZdwdnZucJ2RI0fKut+qqiosXrwY6enpUKlUWLRoEXr27In09HTMmTMH7dq1w4kTJ7BhwwYcOHAAjx49wvDhw/GXv/wFKpUKo0ePxtChQ5GcnIy7d+9i69atGDlyJB49eoSFCxciMzMTgYGBePTokf6asbGxOH36NB49eoSZM2fi7bffBgD4+Phg0aJF+O677/DSSy+hY8eOePnll9GlSxcMGjSoxf/fAgCEjRk8eHCLPn//T0PFhbeHiskfHFQoIiJlZWZmmjsEvTfffFNs2LBBzJ8/X0yaNElUVlYKIYR48sknxU8//SSEECItLU1ERkYKIYQoLCwU1dXVQgghPvnkE7Fq1ao659GZP3++ePbZZ0V1dbXYt2+fcHV1FRcuXBBVVVVi0KBB4ty5c+L27dti5MiR4v79+0IIIdavXy/efvttIYQQ3t7e4sMPPxRCCLFlyxaxePFig9ep7ciRI2LAgAENviIiIhocm56eLp566in966KiIiGEEKNGjRKnT5/Wv19QUKD/fu7cueLrr7/WH6e792+++UZERUUJIYT44x//KBYuXCiEEOL8+fPCwcFBfz7duSorK8WoUaPE+fPn9ff6/vvvCyGEePTokVCr1eKnn34S1dXVYtasWWLSpEkN4jf0b8jY704+URiTcxzwlfcXB1FrePtABjKv3VP0nEGeHfDmlGDJn5s1axYcHBxw//59pKamYtasWfqflZWVAaiZ2vvss8/i+vXrKC8vNzpFc8qUKVCpVNBoNOjevTs0Gg0AIDg4GLm5udBqtcjMzMSIESMAAOXl5YiIiNB/fvr06QCAwYMH4+9//3uT8UdGRuKHH35o1r36+fnhypUrWLFiBSZNmoSxY8caPC45ORkffPABHj58iMLCQgQHB2PKlCkN4svNzQUA/OMf/8DKlSsBACEhIQgJCdGfa8+ePYiPj0dlZSWuX7+OzMxM/c+fffZZADVNUX19fdGnTx8AwNy5c/UNU1uCiYKIFOHi4gIAqK6uRseOHQ3+0l2xYgVWrVqF6OhoHD16FG+99Vaj53vssccAAG3atNF/r3tdWVkJBwcHjBkzptEOD7rPODg4NKtukpycjFdeeaXB++3bt29Q1+jUqRPOnz+PQ4cOYcuWLdizZw/++te/1jmmtLQUy5YtQ3p6Onr27Im33nqrzrTUxuIzNCspJycHf/jDH3D69Gl06tQJCxYsqHMu3X/7xj7fUkwUBmiqszCuIgnANHOHQmSUnL/8Ta1Dhw7w9fXFl19+iVmzZkEIgQsXLmDAgAEoLi6Gl5cXAGD79u36z7i6uuLePWlPRsOGDcPy5cvx73//G/7+/nj48CG0Wi369u3b6GeMXUfKE8WdO3fQtm1bzJgxA71798aCBQv05y8pqVmHpftF3qVLF9y/fx979+7FzJkzjZ73iSeewBdffIHIyEj861//woULFwAA9+7dg4uLC9zc3HDz5k0kJCRg9OjRDT4fEBCAnJwcXL58Gb17926yTVJzcdZTPcUewwEAoypNPwODyFZ98cUX2Lp1KwYMGIDg4GDs378fAPDWW29h1qxZGDlyJLp06aI/fsqUKfjqq68QGhqK48ePN+saXbt2xWeffYbnnnsOISEhGDZsGC5evGj0M3KuY0h+fj5Gjx6N0NBQLFiwAOvWrQNQU9RfunQpQkND8dhjj+E//uM/oNFo8PTTTyM8PLzJ88bGxuL+/fsICQnBBx98gCFDhgAABgwYgIEDByI4OBiLFi3SD7fV5+zsjPj4eEyaNAmPP/44vL29Zd9jbSohDOwOZMVauh9Fxj+/QXXSuwAAzQsbWKMgi5OVlYXAwEBzh0FWzNC/IWO/O/lE0ZQc+X91EBHZAiYKIiIyionCiAyFpx0SEVkjJgoiIjKKiaKeYM8OAGqmyOafTzJzNERE5sdEYcD1LjVTZHvcSeXwExHZPSYKA7wGROHHNpx+SGQI24zXeO+992R/1txdZqViojBCU52FTloOPxHVxjbjNZgoCMcca4af3G5whTZRUxYsWIBVq1YhMjISq1evxuXLlzF+/HgMHjwYI0eO1K+YPnDgAIYOHYqBAwfiqaeews2bN5Gbm4u4uDhs3LhRv2J6wYIFiI2NRWRkJPz8/HDs2DEsWrQIgYGB+nYZQE2b8YiICAwaNAizZs3C/fv3AdS03n7zzTcxaNAgaDQaXLx40eB15FqzZg0ePXqE0NBQzJkzBwDw+eefY8iQIQgNDcWLL76IqqoqVFVVYcGCBejfvz80Gg02btyIvXv36tuRh4aG1mklbrEa7StrpVraZlxc+Yf4V8pBMfmDg+LC20PF/T8NVSYwIoWwzbj524wLIYSLi4v++8zMTDF58mRRXl4uhBAiNjZWbN++vdntyFsb24ybAtuNk6VKWAPc+FHZc3pogAnrJX/MntqM15eUlIQzZ87o+zk9evQI3bp1w5QpU5rVjtzSMVHU5zsSwTgOgLOdiKSwpzbj9QkhMH/+fH1zwNqaakduDZgomiHj2j0EN/6HD5H5yPjL39Tsoc04ADg5OaGiogJOTk6IiorC1KlT8corr6Bbt24oLCxESUkJXFxcmmxHbg1YzG4uNgckajZbbzMOADExMQgJCcGcOXMQFBSEtWvXYuzYsQgJCcGYMWNw/fr1ZrUjt4ZiNtuMG5JzHFP23sN7j96FpjoL1wIXw3PINNYpyCKwzTi1FNuMK4hTZImImCiM+sO9KFxtx7/ciMi+MVE0YlRNvQ13y8wbBxGRubUoUTx48ABVVVVKxWJRJngDms6/vs64dg8nLheYLyCiWmystEitSM6/HUmJorq6Gjt37sSkSZPQrVs3BAQEoEePHggODsZrr72G7OxsyQFYJBatyYI5OzujoKCAyYIkE0KgoKAAzs7Okj4naR1FZGQknnrqKaxbtw79+/dHmzY1eaawsBDJyclYs2YNpk2bhrlz50oKwpI9qAQKS399feJyASJ6d278A0QmplarodVqcfv2bXOHQlbI2dkZarVa0mckJYrvv/8eTk5OyMvL0ycJAHB3d8eMGTMwY8YMVFRUSArAko3yAnCxpk7BYg5ZCicnJ6OtL4iUJun3n5OTEwBg2rRpDX6WlpZW5xhbMMEbcOHadSKyc5ISxZ49e7BmzRqUlJQgKyurTiE7JiZG8sUTExPRr18/+Pv7Y/36hq0IiouLMWXKFP3qzm3btkm+BhERtYykv5dHjBiB0tJSfPrpp1i1ahUuXbqEjh07wtPTE+3atZN04aqqKixfvhyHDx+GWq1GeHg4oqOjERQUpD9my5YtCAoKwoEDB3D79m3069cPc+bMQdu2bSVdi4iI5JOUKLy8vDBv3jz07t1b39q3sLAQOTk5CAgIkHThU6dOwd/fH35+fgCA2bNnY//+/XUShUqlQklJCYQQuH//Ptzd3eHoyLEgIqLWJGnoSTcdT5ckgJpC9uDBg/Uthps7ZS8/Px89e/bUv1ar1cjPz69zzEsvvYSsrCx4enpCo9Fg06ZNdYroOvHx8QgLC0NYWFirzAThegoisieSEkVkZCQ++ugj/Pzzz3XeLy8vx5EjRzB//vw6rYONMZRQVCpVndeHDh1CaGgorl27hh9++AEvvfSSwRbBMTExSE9PR3p6Orp27SrhjhoXPGKS/nvd3tkdbqQpcm4iImsiKVEkJibCwcEBzz33HDw9PREUFARfX1/06dMHu3btwiuvvFJnP1tj1Go1rl69qn+t1Wrh6elZ55ht27Zh+vTpUKlU8Pf3h6+vb5NthJXGxoBEZO8kDfg7Oztj2bJlWLZsGSoqKnDnzh20a9cOHTt2lHzh8PBwZGdnIycnB15eXti9ezd27txZ55hevXohKSkJI0eOxM2bN3Hp0iV9TaO1HHKKQlhpKrx+WXTX4UYa7nkM0w8/cfEdEdk6o08UmZmZdVZZR0VFISMjA0DNeonTp09j8+bNOHXqlOQLOzo6YvPmzRg3bhwCAwPxzDPPIDg4GHFxcYiLiwMAvPHGG0hNTYVGo0FUVBTef//9OpudtAZdc8Cej2qGn4iI7I3RjYt69OiBEydOwMfHBwDQr18/XLp0CQCQmpqKCRMm4Nlnn0VKSgr+53/+x+BCvNamyMZFvzhxuQAdbqQh+VgSXirfigedApEb9gbueQzTH8MnCiKyBbI3Lvruu+/wu9/9Tv+6Q4cO+u937NiBpUuXIj4+HkePHsX777+vULiW55BTFH5sY3hfihOXC/RfRES2yGii0Gg0+OKLL/Sv/f39sXfvXty6dQv79u3D1KlTAQDdunVDWRk3biAiskWSZj1t3LgRf/nLX+Dl5YVBgwZh+PCaGUEVFRW4f/++SQIkIiLzkjTrycPDA4cPH0Z1dXWdhW/JycmIjIxUPDgiIjI/Wf0w6q+OHjt2LMaOHatIQJbqQSVwpRg4lAeM8DB3NERErYfbLDTDKK+aduMPKoFj+Y0fx4I2EdkiJopmmOAN+Ln9ujcFW3kQkT1homhC7TUTAPBjAZCQZ6ZgiIjMQFKNwtXVtUHjPqCmwZ9KpTLYsM+WdHwMwMOa4acRTR5NRGQbJCWKkpISU8VhFdydAc0vC7F1w0/1nziIiGyN7F2AioqKkJ2djdLSUv17TzzxhCJBWYqI3p1ZoCYiuycrUXz66afYtGkTtFotQkNDkZaWhoiICBw5ckTp+IiIyMxkFbM3bdqE06dPw9vbG8nJyTh37pxiGwZZG86AIiJbJytRODs7w9nZGQBQVlaGgIAAfVdZW8PusERk72QlCrVajbt37+Lpp5/GmDFjMHXq1Aa709kil6IsjKtIMjpFljUNIrI1smoUX331FQDgrbfeQmRkJIqLizF+/HhFA7M0xR7D4VKUhakOqfgDonAsv2YhHhGRrZP1RLFx40ZotVoAwKhRoxAdHY22bdsqGpilKVJH4UGnwDpTZImI7IGsRHHv3j2MGzcOI0eOxJYtW3Dz5k2l47J4V4qBNamGh6C4kRER2RJZieLNN99ERkYGtmzZgmvXrmHUqFF46qmnlI7NIrkUZWGFSxL83GqSxbF8znwiItvWol5P3bp1g4eHBzp37oxbt24pFZPFKvao2ahpREUq1g+vaRRIRGTrZCWKjz/+GKNHj0ZUVBTu3LmDTz75BBcuXFA6Noujq1PUphuCeudgBpKy6g7BcQiKiGyBrFlPeXl5+NOf/oTQ0FCFw7Euo7x+/T6v4CGAO4gK7G62eIiITEFWoli/fr3ScVi0ex7DDNYhJnj/OkX21fT2rRwVEVHrkJQoHn/8caSkpDRoN24vbcblqj38xJXeRGRtJCWKlJQUAPbXbjyid2dk3DB+jGN5CSrburZOQERErUj2grv8fCObR9sg7jtBRPZKVo3i3r17GDt2LNzd3TF79mzMnDkT3buziNscHIYiImvDBXcyuBRloZM2qcH7WddLGkyRJSKydlxwJ5Fu0Z3bjdQ67+umyv7z8p3WDomIyKS44E4iQ4vugJppsoE9pBWzuRiPiKyB5BqFEALp6elccEdEZCckP1GoVCqcO3fO7pJERO/OCPbsoH/tXJIHn/R3DdYqiIhsiayhp4iICJw+fVrpWKxGscdwlLp6w7kkr0GtgojI1siaHpucnIy4uDj4+PjAxcVFvzLbnuoUReoo+KS/2+Jz6eoUnCpLRJZKVqJISEhQOg6rpZsqW6SOMncoREQmIStRbN++3eD7v//97yWdJzExEf/5n/+JqqoqLFmyBGvWrGlwzNGjR/Hyyy+joqICXbp0wbFjx+SErJhgzw7IuFbT00q3j7bbjVQmCiKyWbJqFC4uLvovBwcHJCQkIDc3V9I5qqqqsHz5ciQkJCAzMxO7du1CZmZmnWPu3r2LZcuW4euvv0ZGRga+/PJLOeGaTP2pso7lJbIX3XGqLBFZKllPFP/1X/9V5/Wrr76K6OhoSec4deoU/P394efnBwCYPXs29u/fj6CgIP0xO3fuxPTp09GrVy8ANQv8LNkoL+DHgppFd9yXgohsRYtWZus8fPgQV65ckfSZ/Px89OzZU/9arVY3aDT4008/oaioCKNHj8bgwYOxY8cOJcI1GTmL7oiILJ2sJwqNRqPfj6Kqqgq3b9+WXJ8QQjR4r/YeFwBQWVmJM2fOICkpCY8ePUJERASGDRuGvn371jkuPj4e8fHxAIDbt29LikNpjuX21YKdiGyfrERx8ODBX0/g6Iju3bvD0VHaqdRqNa5evap/rdVq4enp2eCYLl266OshTzzxBM6fP98gUcTExCAmJgYAEBYWJvV2Wqz+zKcrxTV7aI/o3YVDUERk9WQNPZ06dQru7u7w9vbGtm3b8Mwzz+Ds2bOSzhEeHo7s7Gzk5OSgvLwcu3fvblDnmDp1Ko4fP47Kyko8fPgQJ0+eRGBgwz5L5lS/SeAoL8DPrWYPbTYIJCJbICtRvPvuu3B1dUVKSgoOHTqE+fPnIzY2VtI5HB0dsXnzZowbNw6BgYF45plnEBwcjLi4OMTFxQEAAgMDMX78eISEhGDIkCFYsmQJ+vfvLydkRdVu5VF/5tMEb2D9cMC7s/Q9tE9cLuDsJyKyOCphqFjQhIEDB+LcuXN4/fXXodFo8Pzzz+vfM7ewsDCkp6eb7gI5xwFAv5YCgH6Fdm7YG/r31qTWDEH5udU8ZYwY2vwd8rhKm4ham7HfnbKeKLy8vPDiiy9iz549mDhxIsrKylBdXd2iIG2NbgjqSjFwzL52jSUiGyMrUezZswfjxo1DYmIiOnbsiMLCQmzYsEHp2CxasGeHOkNQ9emGoPzcWjEoIiITkDXrqX379pg+fbr+dY8ePdCjRw/FgrJoviP1w0/NdaUYeO/vaQCAoYG+Tc6EYqNAIrIkiiy4o8b30dYNQQE1CYMzoYjI2jBRKKCxfbSBX4egOAxFRNZKUqJ44YUXAACbNm0ySTDWKNizAzyHTAO6m3/aLhGRKUiqUZw5cwZ5eXn461//innz5jVow+Hu7q5ocLYor+Ah3jmYAQBcuU1EVkFSoli6dCnGjx+PK1euYPDgwXUShUqlktwY0N6M8gIqb9csxMsreAiAXWaJyPJJGnpauXIlsrKysGjRIly5cgU5OTn6LyaJpk3wBv4QVoLfTw6WtXKbiMgcZE2P/fjjj3H+/HkcP14zTfSJJ55ASEiIooHZA90mR3yqICJLJmvW04cffog5c+bg1q1buHXrFubMmYOPPvpI6dhs2ojeXQBwuiwRWT5ZTxSffvopTp48CRcXFwDA6tWrERERgRUrViganLVxaevQrOM63EhDVOAwJgkisgqyniiEEHBw+PWXooODg8GNiGyW70jFTpV1vQTvHMwwuM82u8kSkSWQ9USxcOFCDB06FNOmTQMA7Nu3D4sXL1Y0MHugG37Kul6zKx5rFURkiWQlilWrVmH06NFISUmBEALbtm3DwIEDlY7N5kUFdkdUYHf9ugoiIkskK1EAwKBBgzBo0CAlY7ErHW6k4Z7Hr3tUGJsBdeJyARsEEpHZsNeTXL4jG9Yqbv7LYGPApnAGFBFZMiYKpfiNBmC4MWBTogK7I7CHq8IBEREpQ1ai2Lx5M4qKipSOxbr1HS+5MWCHG2kmCoaISDmyEsWNGzcQHh6OZ555BomJifY1Nba+esNPziV58El/V9YQlDGcKktE5iIrUaxduxbZ2dlYvHgxPvvsM/Tp0wf//d//jcuXLysdn3XxG41SV284l+TJGoIiIrJEsmsUKpUKHh4e8PDwgKOjI4qKijBz5kz89re/VTI+69J3PFyiN6DU1dvckRARKUbW9NgPP/wQ27dvR5cuXbBkyRJs2LABTk5OqK6uRp8+ffDBBx8oHafV0Q1B6RR7DEeROqrBcb/WKVz1e1VwnwoisiSyEsWdO3fw97//Hd7edf9ybtOmDQ4ePKhIYNZMtzWqjnNJHgAYTBQ6NVNk7zS5T4WuTsF1FUTUWmQNPZWVlTVIEqtXrwYABAYGtjwqK1ekjkJu2Bv6r+YMRUUFduc+FURkkWQlisOHDzd4LyEhocXBWK16M5+CPTsg2LODmYIhIlKWpETx8ccfQ6PR4NKlSwgJCdF/+fr6cuOiJkiZNqurVRjqKKvDqbJE1Fok1Sief/55TJgwAa+//jrWr1+vf9/V1RXu7u6KB2crdDULJWsVREStRVKicHNzg5ubG3bt2mWqeGyKbvgpA1EoUkfVmQXVGHaUJSJLIylRPP7440hJSYGrqytUKpX+fSEEVCoV7t27p3iA9qJ2O4/aXWWJiMxNUqJISUkBAJSUlJgkGCIisjyy96Og5tPPgGrrgAflVeYNhohIIkmJQjfkZKgJoN0PPfmOBHKON3mYS1EWOmmTjBa0iYgsiaREwSGnFvIbDdz8F9xupDJREJHVkLSO4vHHHwdQ82TRoUOHBl/UhF/2rHBp66DIgjy2Hiei1iApUdQuZt+7d6/Bl1SJiYno168f/P3966zLqO/06dNwcHDA3r17JV+DiIhaxmxboVZVVWH58uVISEhAZmYmdu3ahczMTIPHrV69GuPGjTNDlCZy81/AT4mN/lg3VbY5K7SJiExNVqIoLS3F//7v/2L69OmYMWMGNm7ciNLSUknnOHXqFPz9/eHn54e2bdti9uzZ2L9/f4PjPvroI8yYMQPdunWTE6rl+WVvbVw5anT4aUTvLvDu3B55BQ/xz8t3jJ6Sw09EZEqyEsW8efOQkZGBFStW4KWXXkJWVhZeeOEFSefIz89Hz5499a/VajXy8/MbHPPVV19h6dKlcsK0TM3cW5vdZInIUshaR3Hp0iWcP39e/zoyMhIDBgyQdI7GptjW9vLLL+P999+Hg4OD0XPFx8cjPj4eAHD79m1JcZiNbvjpN8ObPpaIyIxkPVEMHDgQaWm/tpw4efIkRowYIekcarUaV69e1b/WarXw9PSsc0x6ejpmz54NHx8f7N27F8uWLcO+ffsanCsmJgbp6elIT09H165dpd2MOdQafmpMhxtp6HAjDY7lJci6XoJ/nkyr0+aDiKi1SHqi0Gg0UKlUqKiowI4dO9CrVy8AwM8//4ygoCBJFw4PD0d2djZycnLg5eWF3bt3Y+fOnXWOycnJ0X+/YMECTJ48GU8//bSk61ikvuP1SULfOPCa4Vljo7yAHwuAY/nABCP7H524XMBd74jIJCQlCiW3OXV0dMTmzZsxbtw4VFVVYdGiRQgODkZcXBwA2FZdojG64ae+4xs9ZIJ3TZIgIjIXSYmi9vanRUVFyM7OrjPbqf72qE2ZOHEiJk6cWOe9xhLEZ599JuncFu+XVdq4ctRootC5UgysSQWGBt7kHhVE1KpkFbM//fRTbNq0CVqtFqGhoUhLS0NERASOHDmidHzWpZn9ngDUGX4CaoagjA0/ATXJAlk5mNYph63IiajVyCpmb9q0CadPn4a3tzeSk5Nx7tw56ygiW6kJ3sD64YCfm/Hj2NKDiExBVqJwdnaGs7MzAKCsrAwBAQG4dOmSooHZjcIcIHGN0ZXaRETmJGvoSa1W4+7du3j66acxZswYdOrUqcHUVrvlO7Lmf5szBKWbJlv4y+yukKbXVLBWQUStTVai+OqrrwAAb731FiIjI1FcXIzx45suyFI9fcfXfCWuadbhhmoVBvWepEx8RESQmShKS0vx5z//GSkpKVCpVHj88cdRXV2tdGxUzwTvmq81qeaOhIjsiaxEMW/ePLi6umLFihUAgF27duGFF17Al19+qWhwVk03BKXT3NlQREQWxmy9nqihplZp1/ZjAZCQZ3y1NhGREmQlCl2vp2HDaubyy+n1RPI12dbD2NNL/ScdIqImmK3Xk91pajbUzX/pi9qd3Ica3VO7qbYe9Z9IlNh2lYjsl9l6PVEtummyAFCYA08AnkOmNWsIiojI1GT3ejp//jyOH6/563jkyJGsUbSEbpos0OypskRErUV2C485c+bg1q1buHXrFubOnYuPPvpI6djIFHKON/wiIjJCVjF769atOHnyJFxcXAAAq1evRkREhH66LFkWY0NYrF8QUVNkPVEIIepsT+rg4GBwa1MyLV07j4Q8c0dCRLZM1hPFwoULMXToUEybNg0AsG/fPixevFjRwGyWlFbkRtRp5wGupyAi05GVKFatWoXRo0cjJSUFQghs27YNAwcOVDo2+6Xb+e43jTcJVLSdR3MSF9dfENktyYlCCAGtVotBgwZh0KBBpojJvtXe+a4Z3WRbylD9gnULIqpNco1CpVLh6aefNkEoBKBmmmz3/pI+omvnQURkCrKK2cOGDcPp06eVjsV+KDiMo6tVGFupTUTUErJqFMnJyYiLi4OPjw9cXFwghIBKpcKFCxeUjs9+FeYg+MJ7AIAH5VX6t4s9htdp79FUOw/F1K5jsF5BZFdkJYqEhASl46Daarf0qMW5pGZ8yVAfKCW7ybJXFBHVJitRdO/evcHGRbGxsUrHZtuMTZOt3dIDQO4vv7h90t81eHiT3WSJiFpAVo1i3rx5yMjIwIoVK/DSSy8hKysLL7zwgtKx0S+a+ot+gjeg6dxKwRCR3eHGRVbGuSQPPunvNqhVEBGZCjcuMieJq7SLPWrWVTRWq9C19BjlZeIhKDkry1kAJ7JashLFyZMnG2xcFBgYqN/YiLOflBfs2QHwnIaMa1EGaxWmbOnBpoJE9k1WokhMTFQ6DpJINwSlEwsgti1wpS1wrGI4AA5LEZEyZCWK2hsYUQs1tUWqAbohKIOnq84DKgEmCiJSiqxEQeZTewjKkOokw1NozU6JDZJY5yAyCyYKapHa9QvWK4hsk2KJ4saNG/Dw8FDqdNQCvtV5UDWyOI/TaolIKlkL7gzhxkWtq7G/3o85Dkem8MaVYui/CktrfuZckge3G0psYEFE9kSxJ4pvvvlGqVNRC5T6ROHd/F+fGK4UA37tgPVhjbcAISIyRtYTxerVq5v1HkngO/LXrxaY4A2sH/7rl5+bQvFZgpzjDb+IyORkPVEcPnwY77//fp33EhISGrxHlqf++gtAuboFd8sjsk2Snig+/vhjaDQaXLx4ESEhIfovX19faDQayRdPTExEv3794O/vj/Xr1zf4+RdffKG/xvDhw+v0l6Lm07X22PZgODKr69YvHIpZtyAi4yQ9UTz//POYMGECXn/99Tq/2F1dXeHu7i7pwlVVVVi+fDkOHz4MtVqN8PBwREdHIygoSH+Mr68vjh07hk6dOiEhIQExMTE4efKkpOvYO11rDwA45BSFQ0516xe7274LPzPERUTWQ1KicHNzg5ubG6ZPnw53d3e4urpi7dq1OHv2LN544w0MHDiw2ec6deoU/P394edX82tq9uzZ2L9/f51EMXz4ryuQhw0bBq1WKyVcQk3NorG+T2tSATxq1XCIyArJqlG8++67mDVrFlJSUnDo0CG8+uqrWLp0qaS/9vPz89GzZ0/9a7VabfTzW7duxYQJEwz+LD4+HvHx8QCA27dvNzsGa6cb/zfWtK85DNUt6rPY9RfGCtpcyU2kCFmznhwcHADUTImNjY3F1KlTUV5eLukcQogG76lUKoPHJicnY+vWrY0Wy2NiYpCeno709HR07dpVUhz27v9VNKxb6NZd6LRk/UXGtXv6LyKyTrKeKLy8vPDiiy/i8OHDWL16NcrKylBdXS3pHGq1GlevXtW/1mq18PT0bHDchQsXsGTJEiQkJKBzZzvZxk3iPhVyjfICjiEKZ2B43YUO118Q2TdZiWLPnj1ITEzEq6++io4dO+L69evYsGGDpHOEh4cjOzsbOTk58PLywu7du7Fz5846x/z888+YPn06/va3v6Fv375yQrULwZ4dZP3Fbqh+saaRB4faw1MWOwxFRCYha+ipXbt2ePDgAXbt2gUAqKioQMeOHSWdw9HREZs3b8a4ceMQGBiIZ555BsHBwYiLi0NcXBwA4J133kFBQQGWLVuG0NBQhIWFNXFWMoVij+Eoda3JKFbVBoQL84gUoRKGigVNiI2NRZs2bXDkyBFkZWWhqKgIY8eOxenTp00RoyRhYWFIT083dxim08gvPaVqAGtSfxl+qrWiu/bWqrqnitywNxS5Xm0mXZzHwjaRUcZ+d8p6ojh58iS2bNkCZ2dnAECnTp0kF7NJWUr9kh3lVTdJXCkGjuUrcmoislKyahROTk6oqqrSz1K6ffs22rRRrBEtmVH9uoWhmkVzptPWx7oGkfWSlShWrlyJadOm4datW/jd736HvXv3Yu3atUrHRobI2Dq1pXQtQABgXMVwjKoGUFzzuuNjgLuz8c87l+QBABMFkZWSlSjmzJmDwYMHIykpCUII7Nu3D4GBgUrHRhIptQCvttotQIC6bUAMTaU1hNNriayb7P0oAgICEBAQoGQsZIGabAGiIFMsytPXbkz9BMZiOdkw7pltrYz9YrpmeZtIyalr1Mc6B5F5MFFQi9SuX9SeRltbscfwhm9KxDoHkfkwUZBstesXV34pbhtKFEXqqBb/gmedg8h8mChs0D2PYY3+rMONNMWuU7t+oXS9gogsBxOFDYroXbd54onLBWaKRFmm3Ma1xSytTQiL66QgJgpSTO16hU5jdQupDNU5mlO3qD2Tivt3E8nDRGFnDA1LKTEcVX+9BWC8biGVoToH6xZErYOJwg7UH4qqTalhKSkty4nIujBRkFUztj7DYuoXRFaOiYKMzpLSkTs81Zx1FnIZW59h9+suLK24TqZl4skLTBR2ztCwlFLDUc1dZyGXsfUZhp4yGmsRwiI3kXFMFGQyXGdBZBuYKKhZag9PKbloj4gsHxMFNWCqBXuG1lnUpnQNo3ahm4VtIvmYKKhVGFpnUZvSNYzahW67L2wTtRATBUnWnFlS9U1AmtEkoHQNo3ahmwvzyObpZrmZaPYTEwU1qTUW7JmTKTZMMhXO0CJzYKIgi9FUDUOOUV5ALJTZOMkY1kDIljFRUKtoarhqaOBNVF6+g0ojxziWl0i6pq7u8Xyvlm+cZAxrIGTrmCioRWoPS7VkGCoqsDuiArsbPUbqtFzd04kSGycZwxoI2bo25g6AiIgsG58oSDGmbAcCyGuRbspeU7WZugai19ahecf5jQb6jjdpKGQ/mCjIZpm615SOseaEZlGYU/O/TBSkECYKslmt1WvK1DWQ2po1PTZxjekDIbvCREEmpVSxuzHGZlOxJxWRMpgoyG40tk7DlLULsynM4ZOFPTFxTYqJglqNsRXegPJPHLWfNkZ5GX66MGXtwmz8Rps7AmpNrVCTYqIgu2BoT2/ARvfJ6DuehWx70gpPjkwUZBcaq2VUts345efBLTo/6yFky5goyGI0NTQlR3OGs/IKHuKdgxktuo5jeYs+3mzPBZXj+aC2rXMxsh66mlT4EiBsoeKnN+vK7MTERPTr1w/+/v5Yv359g58LIbBy5Ur4+/sjJCQEZ8+eNUOUZMtG9O4C787tzR1Gs1wpBvb/u8LcYZCl8RsNuPvWJIsf95rkEmZ7oqiqqsLy5ctx+PBhqNVqhIeHIzo6GkFBQfpjEhISkJ2djezsbJw8eRKxsbE4efKkuUImK9TU9Nzm9JiyFO8czACcnQDfiKYP1u1PQLZPV5MyYa3CbIni1KlT8Pf3h5+fHwBg9uzZ2L9/f51EsX//fsybNw8qlQrDhg3D3bt3cf36dfTo0cNcYZMVa2xoyxb21GjARBvYkIUy8R8GZksU+fn56Nmzp/61Wq1u8LRg6Jj8/HwmClKUKWojptDB2cncIZCl8h0JOHc02enNliiEEA3eU6lUko8BgPj4eMTHxwMALl68iLCwMFkx3b59G127dpX1WWvFe7Y+YZ9I/4y137McdnvPW+T9/svNzW30Z2ZLFGq1GlevXtW/1mq18PT0lHwMAMTExCAmJqbFMYWFhSE9Pb3F57EmvGf7wHu2D6a6Z7PNegoPD0d2djZycnJQXl6O3bt3Izo6us4x0dHR2LFjB4QQSEtLg5ubG4ediIhamdmeKBwdHbF582aMGzcOVVVVWLRoEYKDgxEXFwcAWLp0KSZOnIhvv/0W/v7+aN++PbZt22aucImI7JZZF9xNnDgREydOrPPe0qVL9d+rVCps2bKl1eJRYvjK2vCe7QPv2T6Y6p5VwlDFmIiI6BfcM5uIiIxiovhFU+1ErNXVq1cRGRmJwMBABAcHY9OmTQCAwsJCjBkzBn369MGYMWNQVFSk/8y6devg7++Pfv364dChQ+YKvUWqqqowcOBATJ48GYDt3y8A3L17FzNnzkRAQAACAwNx4sQJm77vjRs3Ijg4GP3798dzzz2H0tJSm7zfRYsWoVu3bujfv7/+PTn3eebMGWg0Gvj7+2PlypUGlx80SpCorKwUfn5+4vLly6KsrEyEhISIjIwMc4eliGvXrokzZ84IIYS4d++e6NOnj8jIyBCvvfaaWLdunRBCiHXr1onf/va3QgghMjIyREhIiCgtLRVXrlwRfn5+orKy0mzxy/XHP/5RPPfcc2LSpElCCGHz9yuEEPPmzROffPKJEEKIsrIyUVRUZLP3rdVqhY+Pj3j48KEQQohZs2aJbdu22eT9Hjt2TJw5c0YEBwfr35Nzn+Hh4SI1NVVUV1eL8ePHi2+//bbZMTBRCCFSU1PF2LFj9a/fe+898d5775kxItOJjo4W3333nejbt6+4du2aEKImmfTt21cI0fDex44dK1JTU80Sq1xXr14VTz75pEhKStInClu+XyGEKC4uFj4+PqK6urrO+7Z631qtVqjValFQUCAqKirEpEmTxKFDh2z2fnNycuokCqn3ee3aNdGvXz/9+zt37hQxMTHNvj6HntB4qxBbk5ubi3PnzmHo0KG4efOmfk1Kjx49cOvWLQC28d/i5ZdfxgcffIA2bX79523L9wsAV65cQdeuXbFw4UIMHDgQS5YswYMHD2z2vr28vPDqq6+iV69e6NGjB9zc3DB27Fibvd/6pN5nfn4+1Gp1g/ebi4kCzW8VYs3u37+PGTNm4E9/+hM6dOjQ6HHW/t/i4MGD6NatGwYPHtys4639fnUqKytx9uxZxMbG4ty5c3BxcTFaa7P2+y4qKsL+/fuRk5ODa9eu4cGDB/j8888bPd7a77e5GrvPlt4/EwWa3yrEWlVUVGDGjBmYM2cOpk+fDgDo3r07rl+/DgC4fv06unXrBsD6/1v885//xNdffw0fHx/Mnj0bR44cwdy5c232fnXUajXUajWGDh0KAJg5cybOnj1rs/f9/fffw9fXF127doWTkxOmT5+O1NRUm73f+qTep1qthlarbfB+czFRoHntRKyVEAKLFy9GYGAgVq1apX8/Ojoa27dvBwBs374dU6dO1b+/e/dulJWVIScnB9nZ2RgyZIhZYpdj3bp10Gq1yM3Nxe7du/Hkk0/i888/t9n71fHw8EDPnj1x6dIlAEBSUhKCgoJs9r579eqFtLQ0PHz4EEIIJCUlITAw0Gbvtz6p99mjRw+4uroiLS0NQgjs2LFD/5lmkVdasT3ffPON6NOnj/Dz8xNr1641dziKOX78uAAgNBqNGDBggBgwYID45ptvxJ07d8STTz4p/P39xZNPPikKCgr0n1m7dq3w8/MTffv2lTQzwtIkJyfri9n2cL/nzp0TgwcPFhqNRkydOlUUFhba9H3//ve/F/369RPBwcFi7ty5orS01Cbvd/bs2cLDw0M4OjoKLy8v8emnn8q6z9OnT4vg4GDh5+cnli9f3mDigzFcmU1EREZx6ImIiIxioiAiIqOYKIiIyCgmCiIiMoqJgoiIjGKiICIio5goiIjIKCYKIonu3r2LP//5z/rXw4cPN8l1tFot/u///s8k5yaSgomCSKL6iSI1NdUk10lKSsLZs2dNcm4iKZgoiCRas2YNLl++jNDQULz22mv4zW9+A6CmjXtAQACWLFmC/v37Y86cOfj+++8xYsQI9OnTB6dOndKf4/PPP8eQIUMQGhqKF198EVVVVXWukZKSglWrVmHv3r0IDQ1FTk5Oq94jUR0KtiQhsgv1N5FxcXHRv+/g4CAuXLggqqqqxKBBg8TChQtFdXW12Ldvn5g6daoQQojMzEwxefJkUV5eLoQQIjY2Vmzfvr3BdcaNGyd+/PFH098QURMczZ2oiGyJr68vNBoNACA4OBhRUVFQqVTQaDTIzc0FUDOkdObMGYSHhwMAHj16pG8TXdulS5fQr1+/VoudqDFMFEQKeuyxx/Tft2nTRv+6TZs2qKysBFDT+n3+/PlYt25do+cpKCiAm5sbnJycTBswUTOwRkEkkaurK0pKSmR/PioqCnv37tVvX1lYWIi8vLw6x+Tk5Fj1xjpkW5goiCTq3LkzRowYgf79++O1116T/PmgoCCsXbsWY8eORUhICMaMGaPfrUwnICAAd+7cQf/+/U02q4qoubgfBRERGcUnCiIiMoqJgoiIjGKiICIio5goiIjIKCYKIiIyiomCiIiMYqIgIiKjmCiIiMio/w92pSbw4fLAYwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for treatment_type in (\"standard\", \"test\"):\n",
    "    mask_treat = data_x[\"Treatment\"] == treatment_type\n",
    "    time_treatment, survival_prob_treatment, conf_int = kaplan_meier_estimator(\n",
    "        data_y[\"Status\"][mask_treat],\n",
    "        data_y[\"Survival_in_days\"][mask_treat],\n",
    "        conf_type=\"log-log\",\n",
    "    )\n",
    "\n",
    "    plt.step(time_treatment, survival_prob_treatment, where=\"post\", label=f\"Treatment = {treatment_type}\")\n",
    "    plt.fill_between(time_treatment, conf_int[0], conf_int[1], alpha=0.25, step=\"post\")\n",
    "\n",
    "plt.ylim(0, 1)\n",
    "plt.ylabel(r\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the results are inconclusive, because the difference between the two estimated survival functions is too small to confidently argue that the drug affects survival or not.\n",
    "\n",
    "*Sidenote: Visually comparing estimated survival curves in order to assess whether there is a difference in survival between groups is usually not recommended, because it is highly subjective. Statistical tests such as the* [log-rank test](https://scikit-survival.readthedocs.io/en/latest/api/generated/sksurv.compare.compare_survival.html#sksurv.compare.compare_survival) *are usually more appropriate.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Survival functions by cell type\n",
    "\n",
    "Next, let's have a look at the cell type, which has been recorded as well, and repeat the analysis from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x120eb8190>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAABEwElEQVR4nO3deVhTd/Y/8HcCKIuAAiJItIAg+04EKipoBUWLC9raugxai6hdrVbbzlTbsdV2Zro4boO1Wr91a+0o/bnQurbiggoqVkQpgqwqi+yCBO7vj0yuhISQXBJCwnk9D89Dbm7uPZdaDvd+Pud8eAzDMCCEEEI6wNd2AIQQQno2ShSEEEIUokRBCCFEIUoUhBBCFKJEQQghRCFKFIQQQhTSaqJYsGABbG1t4e3tLfd9hmHwxhtvwMXFBb6+vsjIyOjmCAkhhGg1UcTHxyMlJaXD948dO4acnBzk5OQgKSkJixcv7sboCCGEAFpOFKNHj4aVlVWH7ycnJ2PevHng8XgIDQ1FVVUVSktLuzFCQgghPXqMori4GEOGDGFfCwQCFBcXazEiQgjpfQy1HYAi8rqL8Hg8mW1JSUlISkoCAGRnZ8Pd3V3jsRFCiD7Jz89HeXm53Pd6dKIQCAQoLCxkXxcVFWHw4MEy+yUkJCAhIQEAEBwcjCtXrnA+5+Wr2yEMeAUAUJ92CWYhIzgfixBCdEVwcHCH7/XoR0+xsbHYtWsXGIbBxYsXYWlpCXt7e22HRQghvYpW7yheeuklnDlzBuXl5RAIBPjoo4/Q3NwMAEhMTERMTAyOHj0KFxcXmJqaYseOHRqPKa+8HqLcCoQNs9b4uQghRBdoNVHs3btX4fs8Hg+bNm3qpmgAHFuFiFvH8GdtETDscwD0+IkQQnr0GIU22Nbko7U4VdthEKJWzc3NKCoqQmNjo7ZDIVpmbGwMgUAAIyMjpT9DiaKtievx8NYxbUdBiNoVFRXB3Nwcjo6OcmcOkt6BYRhUVFSgqKgITk5OSn+uRw9md7fPLn2GZRatOMqvw/1PP8WDdetQd/q0tsMipMsaGxthbW1NSaKX4/F4sLa2VvnOku4o2rluBAAMYgA0ZWcDEI9TdITGL4iuoCRBAG7/DuiOoo2VI1bCTzzpCnbvv4++VLhHCFHCV199hV27dmns+H/729/g6+sLf39/REVFoaSkBIC4SM7ExAT+/v7w9/dHYmIi+5nnnnsOjx49Usv5KVF04EJuhbZDIIToAJFIhG+//RYvv/yyxs6xYsUKZGZm4tq1a5g8eTI+/vhj9r1hw4bh2rVruHbtGrZu3cpunzt3LjZv3qyW81OiIIRoXH19PSZNmgQ/Pz94e3tj//79AICUlBS4u7sjPDwcb7zxBiZPngwAWLNmDf75z3+yn/f29kZ+fj4AYOrUqQgKCoKXlxfbugcA+vXrh5UrVyIoKAjPPfccLl26hIiICDg7O+Pnn38GIB6rmT9/Pnx8fBAQEIDT/xuD3LlzJ1577TX2WJMnT8aZM2fQ0tKC+Ph4eHt7w8fHB19++aXMtZ06dQqBgYEwNBQ/yY+IiMDKlSsxYsQIDB8+HGfPnu3yz8/CwkLqZ6nM46PY2NhOSxCURWMUbXz0/26iqUWcO0/eegDN/X1ASO+SkpKCwYMH48iRIwCA6upqNDY24tVXX8WpU6fg4uKCF198Ualjffvtt7CyssLjx48hFAoRFxcHa2tr1NfXIyIiAp999hmmTZuGv/71rzh+/DiysrLwl7/8BbGxsWxd1o0bN5CdnY2oqCjcuXOnw3Ndu3YNxcXF+OOPPwAAVVVVMvucO3cOQUFBUttEIhEuXbqEo0eP4qOPPsKJEyek3q+trcWoUaPknnPPnj3w9PSU2f7BBx9g165dsLS0ZBMcAOTl5SEgIAAWFhZYu3Yte9wBAwagqakJFRUVsLbuWgExJQo5GAY49UeWUolCMtBNg9pEV3z0/24iq6RGrcf0HGyB1c97dfi+j48Pli9fjpUrV2Ly5MkYNWoUrl27BicnJ7i6ugIA5syZI3WH0JENGzbg4MGDAIDCwkLk5OTA2toaffr0wYQJE9jz9e3bF0ZGRvDx8WHvRlJTU/H6668DANzd3fHMM88oTBTOzs64e/cuXn/9dUyaNAlRUVEy+5SWlsLDw0Nq2/Tp0wEAQUFB7LnbMjc3x7Vr1zq91rY++eQTfPLJJ1i3bh02btyIjz76CPb29igoKIC1tTXS09MxdepU3Lx5k70DsbW1RUlJSZcTBT16amP1817oa9AKyV1dQ5NIuwERoieGDx+O9PR0+Pj44L333mOfsXf0CMXQ0BCtra3sa8l0zjNnzuDEiRO4cOECrl+/joCAAPY9IyMj9nh8Ph99+/ZlvxeJxP8vy+tIreh8AwYMwPXr1xEREYFNmzZh4cKFMp81MTGRmW4qObeBgQF77rZqa2vZAej2X1lZWXJjlHj55Zfx008/seeRJIGgoCAMGzZMKvE1NjbCxMRE4fGUQXcUcmQbt6KfSToA8RTZutOn0S8yUstREaIeiv7y15SSkhJYWVlhzpw56NevH3bu3Il3330XeXl5yM3NxbBhw6Sepzs6OuLw4cMAgIyMDOTl5QEQP7IaMGAATE1NkZ2djYsXL6oUx+jRo7F7926MHTsWd+7cQUFBAdzc3FBTU4PNmzejtbUVxcXFuHRJ/KSgvLwcffr0QVxcHIYNG4b4+HiZY3p4eODPP/9UKQ5V7yhycnLYO6+ff/6ZXUqhrKwMVlZWMDAwwN27d5GTkwNnZ2cA4qR4//59ODo6qhSbPJQo2hnbxMN1IwZNJjfQ7B8Mg7w/UX/xIiUKQrrgxo0bWLFiBfh8PoyMjLBlyxYYGxsjKSkJkyZNgo2NDcLDw9mxgLi4OOzatQv+/v4QCoUYPnw4AGDChAnYunUrfH194ebmhtDQUJXiWLJkCRITE+Hj4wNDQ0Ps3LkTffv2xciRI+Hk5AQfHx94e3sjMDAQgHjxtPnz57N3G+vWrZM55sSJEzF37tyu/Hg6tWrVKty+fRt8Ph/PPPMMO7vp999/x4cffghDQ0MYGBhg69at7Kqh6enpCA0NZQfZu4LHdHQvpqO6uh7Fwy/8sNQMyOUNwX8Gvgnr78TTywa9957Cz9EYBenJbt26JfMcvac5c+YM/vnPf7J3Erpk2rRp+Pzzz9m/+nuCN998E7GxsRg3bpzMe/L+PSj63UljFHL0YxpgxVShoO5PVDbUoqFJhKySaoWfUVS9TQjRb+vXr0dpaam2w5Di7e0tN0lwQY+e2ikY6AvUX8QA1MK4Nl/b4RDSa0RERCAiIkLbYXDi5uYGNzc3bYch5dVXX1XbseiOop27dkLU8Uy1HQYhhPQYlCgIIYQoRIlCBVkl1Z2OVRBCiL6hRNEJfgutCEYI6d1oMLsDrQyQWmeHmQD4pcUwSdrAvveg79Mfm1loKFtjQetrE0L0Ed1RyGHAE5eWXKm3RZ2bA1rtHeTu96SgAPUqVoYSQtSnX79+AMTrMnh7e3M6Rnx8PA4cOABAPPOqo1qCGTNm4O7du9wCVUFsbKzUtRQUFCAyMhIBAQHw9fXF0aNHAYirsiW9rTSN7ijkMOQzaGHEPWNqvR1hMv45qfedBlsCAB7IqdIkhOifmzdvoqWlhW2PoSn//e9/2eQnsXbtWrzwwgtYvHgxsrKyEBMTg/z8fAwcOBD29vY4d+4cRo4cqdG46I6iA/2YBsQgVdthEKIXOlqPwtHREe+//z7CwsIQHByMjIwMREdHY9iwYWybirq6OowbNw6BgYHw8fFBcnKywnO1tLRg+fLl8PHxga+vL/79738DELe0GDNmDIKCghAdHa1Sgdzu3bsxZcoU9nW/fv3wwQcfwM/PD6GhoXjw4IGqPxIZdXV1+OKLL/DXv/5VajuPx0NNjbjbb3V1NQYPHsy+N3XqVOzevbvL5+4M3VHI0WhkBiNRI8biMoqgnspGQnqMY6uA+zfUe0w7H2Di+g7flrcehcSQIUNw4cIFvP3224iPj8e5c+fQ2NgILy8vJCYmwtjYGAcPHoSFhQXKy8sRGhqK2NjYDjvPJiUlIS8vD1evXoWhoSEqKyvR3NyM119/HcnJyRg4cCD279+PDz74AN9++61Sl3fu3Dm89NJL7Ov6+nqEhobik08+wbvvvott27bJ/II/ffo03n77bZljmZqa4vz58zLb//a3v+Gdd96Bqal0HdeaNWsQFRWFf//736ivr5da2yI4OFjmvJpAiUKOxj4W4DU+1nYYhOgNeetRSMTGxrL71NXVwdzcHObm5jA2NkZVVRXMzMzw/vvv4/fffwefz0dxcTEePHgAOzs7uec6ceIEEhMT2WZ4VlZW+OOPP/DHH39g/PjxAMR3Hfb29krHX1paioEDB7Kv+/Tpw67GFxQUhOPHj8t8JjIyUukOsdeuXcOff/6JL7/8Umb9ir179yI+Ph7vvPMOLly4gLlz5+KPP/4An89n15vQNEoUSiioE7cQHtrPpdN9ufZ8otlSpNso+MtfUyTrURw9ehTvvfceoqKi8OGHHwKA1LoRku8lr0UiEXbv3o2ysjKkp6fDyMgIjo6OMus/tMUwjMzdBsMw8PLywoULFzjF337NibZrX3S05oQqdxQXLlxAeno6HB0dIRKJ8PDhQ0RERODMmTPYvn07UlJSAABhYWFobGxEeXk5bG1t1bbeRGdojKKdBqunSxA+ZgyRWif/rxZCiPJKSkpgamqKOXPmYPny5cjIyFD6s9XV1bC1tYWRkRFOnz6Ne/fuKdw/KioKW7duZX95V1ZWws3NDWVlZWyiaG5uxs2bN5WOgcuaE5I7ivZf8h47LV68GCUlJcjPz0dqaiqGDx+OM2fOAACGDh2KkydPAhB3fW1sbGTvbu7cucN5tpcqKFF0INu4Fb+bN+NKva3Me5IK7YYmEZ4UFODBunWoa7OGLSFE2o0bNzBixAj4+/vjk08+Uem5+uzZs3HlyhUEBwdj9+7d7KI9HVm4cCGGDh0KX19f+Pn5Yc+ePejTpw8OHDiAlStXws/PD/7+/nJ/YXdk0qRJ7C/u7vavf/0L27Ztg5+fH1566SXs3LmTvZs5ffo0Jk2apPEYaD2Kdnamn0B24U84Up0K98c8VJWuRKKLeMZB+0dPhpfOwfLWdTwpKECfoUM7XbNCEXr0RDRJF9aj6MkeP36MyMhInDt3DgYGBtoOhzV69GgkJydjwIABKn2O1qNQg0AzD/g1d76faMRIDHrvPfQZOlTzQRFCtMbExAQfffQRiouLtR0Kq6ysDMuWLVM5SXBBg9mEEKKE6OhobYcgZeDAgZg6dWq3nKtLdxT19fVoaWlRVyw9Tj88liq6K6j7k50BpW60Qh4hpKdSKVG0trZiz549mDRpEmxtbeHu7g57e3t4eXlhxYoVyMnJ0VSc3a7RyAwAMJLJwNY/LbA5zxYXKs1k9qO244QQfadSooiMjERubi7WrVuH+/fvo7CwEA8fPsTZs2cRGhqKVatW4fvvv9dUrN2qsY8FGnnG6MsT3zEVN/bB1WrZREEIIfpOpTGKEydOwMjICPfu3QOf/zTHWFlZIS4uDnFxcWhuVmIUWEcY8Foh6FOHtwZk4otympVECOmdVLqjMDIyAgBMmzZN5r2L/2u3LdlHH/xpxMObAw2QZljQ6b5N2dlUS0FIF+3cuROvvfaaRs/BMAzGjh3LNtpTt8LCQkRGRsLDwwNeXl74+uuv2fdefPFF+Pv7w9/fH46OjvD39wcgrjOJj4/XSDzqoFKi+OGHH7Bq1SrU1tbi1q1bUgPZCQkJKp88JSUFbm5ucHFxwfr1sm0Fqqur8fzzz8PPzw9eXl7YsWOHyufgyttkGFyaGdw1BLL5isdezEJDAYDWpiBEBxw9ehR+fn6wsLDQyPENDQ3xr3/9C7du3cLFixexadMmZGVlAQD279/PVmjHxcVh+vTpAMR9roqKilBQ0PkfpdqgUqIYOXIkPD098ejRIyxbtgyurq4IDAzE5MmTVe430tLSgqVLl+LYsWPIysrC3r172R+mxKZNm+Dp6Ynr16/jzJkzeOedd/DkyROVzsNVoJkH3q11htuTJ+jTKt0gsP3MpwK3QLQ4uaChScRWbbf9IoSIW2IHBQXBy8sLSUlJ7PYdO3Zg+PDhGDNmDM6dO8duLysrQ1xcHIRCIYRCIfvemjVrsGDBAkRERMDZ2RkbNjxdffKLL76At7c3vL298dVXX8mNo23L8Pz8fHh4eODVV1+Fl5cXoqKi8Phx1xqC2tvbIzAwEABgbm4ODw8PmfoLhmHwww8/SHWkff7557Fv374unVtTVBqjcHBwwLx58zBs2DB2oYzKykrk5eV1Wlbf3qVLl+Di4sIuBDJr1iwkJyfD0/NpryUej4fa2lowDIO6ujpYWVmxHSG7Q6GJJ5r5T1e0Km7sg815tgiwrMfQfgo+SEgP9tmlz5Bdma3WY7pbuWPliJUK9/n2229hZWWFx48fQygUIi4uDk+ePMHq1auRnp4OS0tLdiU3AHjzzTfx9ttvIzw8HAUFBYiOjsatW7cAANnZ2Th9+jRqa2vh5uaGxYsXIzMzEzt27EBaWhoYhkFISAjGjBnDHk/i3Llz+M9//sO+zsnJwd69e7Ft2za88MIL+OmnnzBnzhypz+zevRv/+Mc/ZK7JxcWFXR1Pnvz8fFy9ehUhISFS28+ePYtBgwbB1dWV3RYcHIz169fj3XffVfhz1AaVfutKujK2XU3JysoKVlZWMvt0pri4GEOGDGFfCwQCpKWlSe3z2muvITY2FoMHD0ZtbS32798vNYgukZSUxP6FUlZWpsolKS3Asl4cd2MfAMCLGjkLIfprw4YNOHjwIADxc/ycnBzcv38fERERbJO7F198EXfu3AEgnjzT9ilDTU0NamtrAYh7L/Xt2xd9+/aFra0tHjx4gNTUVEybNg1mZuLZidOnT8fZs2dlEkVlZSXMzc3Z105OTuxYQVBQkEybb0Dcb2r27NkqXW9dXR3i4uLw1VdfyTzm2rt3r9TdBIBuaxnOhUqJIjIyEnFxcZgyZQqGtmlb8eTJE6SmpuK7775DZGSkUoMy8lpMtU8wv/zyC/z9/XHq1Cnk5uZi/PjxGDVqlMwPPSEhgR0jCQ4OVuWSlBZmVY8wq3pszpNtEkiILunsL39NOHPmDE6cOIELFy7A1NQUERERbNvujv6wbG1txYULF+Q+1m7bjlzS5lvZtnWGhoZobW1l/+hsfyx5j55UvaNobm5GXFwcZs+ezY5DSIhEIvz3v/9Fenq61PbuahnOhUpjFCkpKTAwMMBLL72EwYMHw9PTE05OTnB1dcXevXvZFaqUIRAIUFhYyL4uKiqSWuIPED+7nD59Ong8HlxcXODk5ITsbPXeMsvTttU4IaTrqqurMWDAAJiamiI7O5udJRkSEoIzZ86goqICzc3N+PHHH9nPREVFYePGjezrzhYBGj16NA4dOoSGhgbU19fj4MGDUgskSbi5ueHu3btyjtCx2bNny20ZLi9JMAyDV155BR4eHli2bJnM+ydOnIC7uzsEAoHU9u5qGc6FSncUxsbGWLJkCZYsWYLm5maUl5fDxMQE/fv3V/nEQqEQOTk5yMvLg4ODA/bt24c9e/ZI7SPpwz5q1Cg8ePAAt2/f1vji5l4OlgCAvEqNnoaQXmXChAnYunUrfH194ebmhtD/zRS0t7fHmjVrEBYWxg4CS2ZTbtiwAUuXLoWvry9EIhFGjx7NrqMtT2BgIOLj4zFihLjmaeHChTKPnYCnLcNdXDpfiIyLc+fO4f/+7//g4+PDPtL69NNPERMTAwDYt2+fzGMnoPtahnOhsM14VlYWPv30U7baety4cdiwYQO8vLwAAD///DMyMzMRFRXF/sdRxdGjR/HWW2+hpaUFCxYswAcffMD+Q0hMTERJSQni4+NRWloKhmGwatUqmUGm9rraZvzy/csAgLwbF2CWU4y9hocBANNMxHO7JY+e1vtIP/4ySdoAg7w/0TjtRYhGjERHPAdbdvgetRonmkJtxp8qLS3FvHnz5C5fqi1NTU0YM2YMUlNTu2XCjqptxhVGNG7cOKmlA4uKitgkcf78ecydOxcvvvgi4uPj8cknn8gtxFMkJiaGzbISiYmJ7PeDBw/Gr7/+qtIxNcW4Nh+N5o4dvt/sHwSDvD9hdC1dYaJQRNIYkBIGIZpjb2+PV199FTU1NRqrpVBVQUEB1q9f362zOlWhcIzi119/xQcffMC+bvtD3bVrFxITE5GUlIQzZ87gs88+01yUOkA0YiRanDRzK0sIUa8XXnihxyQJAHB1dUVERIS2w+iQwkTh4+OD3bt3s68lI/wPHz7EoUOH2KIVW1tbNDU1aTZSQgghWqHSrKcvv/wS//nPf+Dg4IDAwEA8++yzAMRTwerq6jQSYHcT2gnZ7+tdHbQYCSGE9AwqPRCzs7PD8ePHpeYgA+LR+sjISLUHp4+Uaekh7HQPQgjpPpxWuGtfHR0VFSXVu4UQQtrr16/7+94cOnQIH3/8scaO/8UXX8DT0xO+vr4YN24c7t27B0D8x7OkS6y/vz+MjY1x6NAhAOJ2Rbq2yFuXlkIlhBBNYBgGra2tXT7O559/jiVLlqghIvkCAgJw5coVZGZmYsaMGWyfpsjISLYo79SpUzA1NUVUVBQAYPHixfj88881FpMm9My5WDqgbQfZof3UO9tJ0frZNHWW6Lq6ujpMmTIFjx49QnNzM9auXYspU6YgPz8fEydORGRkJC5cuIBDhw5h165d2L17N4YMGQIbGxsEBQVh+fLlyM3NxdKlS1FWVgZTU1Ns27ZNpjHpnTt30LdvX9jY2AAA4uPjYWFhgStXruD+/fv4/PPPMWPGjC5dS9tH7qGhoXJX+Dxw4AAmTpwIU1NTAMCoUaMQHx8PkUjUY6fDtqcbUWpBg5UnTCvFDcn+NOLhP4YX4cXUA7Blu8gCkOkkyy8thknSBjlHVKzZP4hz/QUhusTY2BgHDx6EhYUFysvLERoaitjYWADA7du3sWPHDmzevBlXrlzBTz/9hKtXr0IkEiEwMBBBQUEAxP3dtm7dCldXV6SlpWHJkiU4deqU1HnOnTvHtvuWKC0tRWpqKrKzsxEbGys3UYwaNYptPtjWP//5Tzz33HMdXtf27dsxceJEme379u2TauXB5/Ph4uKC69evs9fT06mUKMzNzeU28JJ0jNXUilHa9FxDK2DKx+0+NWhlchBg+fQv+vadZJv9g8BlfT9+aTGMAEoUpFvc//RTNN1Sb8+0vh7usHv/faX2ZRgG77//Pn7//Xfw+XwUFxfjwYMHAIBnnnmGbe+RmpqKKVOmsI3ynn/+eQDiO5Lz589j5syZ7DHlTc8vLS1lu9JKTJ06FXw+H56enuw52zt79qxS19HW999/jytXruC3336TieHGjRuIjo6W2i7pFKuXiUJeltVHQjshbhafAADEtPZBfEk+XhC4osHoaRdZADKdZEUjRnL6Zc/lDoQQXbV7926UlZUhPT0dRkZGcHR0ZDvJSlqEA/I7TAPirrL9+/fvtEmgiYkJqqulZxm27RTb0fFVvaM4ceIEPvnkE/z2229SxwfEq4JOmzZNZonontwpVh7Oj54ePXqEnJwc9j8wIO7eqG8KBvrCtiYffVofo0HbwRCiBsr+5a8p1dXVsLW1hZGREU6fPs3OFGovPDwcixYtwnvvvQeRSIQjR47g1VdfhYWFBZycnPDjjz9i5syZYBgGmZmZ8PPzk/q8h4eH3DGDzqhyR3H16lUsWrQIKSkpsLWVXYJg7969WLduncz2O3fusO2QdAGnWU/ffPMNRo8ejejoaKxevRrR0dFYs2aNmkPTLkkX2bt2QlQY2QMA+C2Nij5CCFHC7NmzceXKFQQHB2P37t0dro4pFAoRGxsLPz8/TJ8+HcHBwbC0FP9/uXv3bmzfvh1+fn7w8vJCcnKyzOdHjx6Nq1evKr1OBRcrVqxAXV0dZs6cCX9/f3asBRCvbldYWIgxY8ZIfebBgwcwMTGBvb29xuJSN053FF9//TUuX76M0NBQnD59GtnZ2Vi9erW6YyOE6BFJ9wYbGxupZqNt/fHHH1Kvly9fjjVr1qChoQGjR4/GO++8A0C8Kl1KSorC85mamuK5557DyZMn8dxzz2Hnzp1y4+mKEydOdPieo6OjzFrZALBnzx4sWrSoy+fuTpzuKIyNjWFsbAxAPIjk7u6O27dvqzUwXZHbYIz9BaXddr76tEvsFyH6LiEhAf7+/ggMDERcXJzMLKbOvP/++2ho6FkPjfv374+//OUv2g5DJZzuKAQCAaqqqjB16lSMHz8eAwYMkFmdrjcIsKxHboMxrlab0RrahGhA+8XMVDVo0CCpx0E9wfz587Udgso4JQrJAulr1qxBZGQkqqurMWHCBLUG1hOV8muwo/UQfHiuCOZ5IcyqHlerzTr/oBIM8v6E4aVzEI0YiaySaoULHBFCSHfi9Ojpyy+/RFFREQBgzJgxiI2NRZ8+fdQaWE/zXEMr7FstcB/luMGot09Ls794LrXRtfRO9iSEO00O6hLdweXfAadEUVNTg+joaIwaNQqbNm3qsHBFnzxfz2BRYyjsYKP2Y9OiR0TTjI2NUVFRQcmil2MYBhUVFewYs7I4PXpavXo1Vq9ejczMTOzfvx9jxoyBQCBQOANAF0naeLSYiIto6l0d0FIlHrSvdXYV75TXedtwQrRNIBCgqKgIZWVl2g6FaJmxsTEEAoFKn+lSrydbW1vY2dnB2toaDx8+7MqheiQvB0vkVUpvMxA1oMXQVDsBEcKRkZERnJyctB0G0VGcEsWWLVuwf/9+lJWVYcaMGdi2bRs8PT3VHZvOKG7sg1U3atDXQH5JfrgNH+NtDdR+3vZTZKmzLCFEEzglinv37uGrr76Cv7+/msPRPQGW9Qrfz29ggPJWjSQKQgjpDpwSxfr169Udh85q2yRQ3roUq7OaOR2XpsgSQnoKlWY9hYeHAxC3G7ewsGC/JK/1UYOV+JGabU0+nO9fBgAUNeQis/K8NsMihJBuo1KiSE1NBSBuN15TU8N+SV7rq4KBvgCAoWWZ8DYZBgC4VZWhzZAIIaTbcC64k9fsSl/dtRPioYUjACDQzAMC02HaDYgQQroRpzGKmpoaREVFwcrKCrNmzcKMGTMwaNAgdcfWI3g5WAJ8MxgbGaCxuYXTMfIbGHasQlMzoAghRFM43VGsXr0aN2/exKZNm1BSUoIxY8YoXEtWFwnthErtV+vpilpP1w7fD7fhw9FUvHxsfgOD1PJWtcQnD3WVJYRoAhXccWAgErctNq3MYrcZ1xYDcmY9jbc1YO8gVJ0BlVXytOqbZkARQrSF0x3Fli1bEBERgXHjxqG8vBzbtm1DZmamumMjhBDSA6h8R8EwDK5cuUIFd4QQ0kuofEfB4/Fw9epVShIawC8thknSBhheOqftUAghhMXp0VNYWBguX76s7lh6tWb/ILTaO4BfWkzrUhBCehROg9mnT5/G1q1b4ejoCDMzMzAMAx6Pp/fjFMZGBnCyMQPKuR+j7VRZKf1GAOEjsPjXzbBsZiC/vSAhhHQ/Toni2LFj6o6jVwi34QOdTI9tbGEA8ChREEJ6DE6J4rvvvpO7/cMPP1TpOCkpKXjzzTfR0tKChQsXYtWqVTL7nDlzBm+99Raam5thY2OD3377jUvI3aKg7s8O3xvaz0VqqmxH7v/Kk7u97VRZCZoySwjpDpwShZmZGft9Y2MjDh8+DA8PD5WO0dLSgqVLl+L48eMQCAQQCoWIjY2VWteiqqoKS5YsQUpKCoYOHdrttRqSorvLpZ0/Uqt3dYDBfU1HRAgh3Y9TonjnnXekXi9fvhyxsbEqHePSpUtwcXGBs7MzAGDWrFlITk6WShR79uzB9OnTMXToUADiAj9tsqjIxcD8C1qNgRBCuhunWU/tNTQ04O7duyp9pri4GEOGDGFfCwQCmUaDd+7cwaNHjxAREYGgoCDs2rVLHeFyUikIBABYFVHXWEJI78LpjsLHxwc8nvhZektLC8rKylQen2AYRmab5JgSIpEI6enpOHnyJB4/foywsDCEhoZi+PDhUvslJSUhKSkJADS2eHyZY5hUkrjfXIFd5YcBAN4mwxBoptqjN0UcKovRJ2mD3Pea/YMgGjFSbecihJDOcEoUhw8ffnoAQ0MMGjQIhoaqHUogEKCwsJB9XVRUhMGDB8vsY2NjAzMzM5iZmWH06NG4fv26TKJISEhAQkICACA4OFjVy1GZZE0KQJwwAKgtUVx1DAAAOMl5j19aDCOAEgUhpFtxShSXLl3ChAkTYG5ujrVr1yIjIwN//etfERgYqPQxhEIhcnJykJeXBwcHB+zbtw979uyR2mfKlCl47bXXIBKJ8OTJE6SlpeHtt9/mErJaBZp5sIlBclcBAMa1+Wg0d+zSsS8OD8PF4WH4yNNI5j2TDu4yCCFEkzglir///e+YOXMmUlNT8csvv2D58uVYvHgx0tLSlD+xoSE2btyI6OhotLS0YMGCBfDy8sLWrVsBAImJifDw8MCECRPg6+sLPp+PhQsXwtvbm0vIeknelFkAEOVWqHScsGHW6giHEKKnOCUKAwNxLcCRI0ewePFiTJkyBWvWrFH5ODExMYiJiZHalpiYKPV6xYoVWLFiBZcw1cfeF1BiiqwikhqLoXJakaubYab0gLvIV/k7PUIIaY/TrCcHBwcsWrQIP/zwA2JiYtDU1ITWVs0tyEMIIUR7OCWKH374AdHR0UhJSUH//v1RWVmJf/zjH+qOTefUuzpoOwRCCFE7To+eTE1NMX36dPa1vb097O3t1RZUT2ZRkQvn+5dxV8mlUgkhRNeppeCut5AU3Q2vugknGzNxJ1kNyKplcPxhi0aOTQghqqJEoYIyxzDUWA+T2X7vSSky6m+p5RzhNuL/JKkddJmVLG5ECxwRQrqLSoli7ty5AICvv/5aI8HoIknx3R+Pc9VyvPG2BvA0l99BVrK4EQC1LnB0IbdC4RchpHdTaYwiPT0d9+7dw7fffot58+bJtOGwsrJSa3C6INDMQ21JojOiESPZqmwqviOEdBeVEkViYiImTJiAu3fvIigoSCpR8Hg8lRsD9kYFdX92Sy1FW+3rKtqiGgtCSGdUevT0xhtv4NatW1iwYAHu3r2LvLw89ouSBCGE6CdO02O3bNmC69ev4+zZswCA0aNHw9fXV62BEUII6Rk4JYoNGzYgKSmJraWYPXs2EhIS8Prrr6s1uJ5CaCcEHjfictUdmfecbMyA8qevW+ya0WDlKbNfW6L+ATLbFD0eIoQQbeKUKL755hukpaWxS6KuXLkSYWFhepso2jOtLoZb6ib2dX+TZgCAs0h9hXiSWorO1tgmhBBN45QoGIZhGwMC4iaB8hYi0keSorv2rhsBmTXX0U8NiSLcho+s2haklrcqTBSSmgoJTS1q1H6KLHWbJaR34ZQo5s+fj5CQEEybNg0AcOjQIbzyyitqDaynKnMMQ5ljmNQ2wb0ruF6dilN9Gai2crh8420NOiy4k2j2D0LbFStoUSNCiKZwShTLli1DREQEUlNTwTAMduzYgYAA2efuesVpFHBVdowCENdSFJWndms4bWsqAKqrIIRoDqdEAQCBgYEqrWjXm5hWZnU6oE0IIbqCc6Ig3N2uugoAcJMz+6m7yZttRUV4hJC2KFGogZONGfgFPPZ7AIC9Jfv+zWL5S5YSQogu4NQ9duPGjXj06JG6Y+nxhP2HQ9h/uEaOLe+v+PwGBquzmrE6q5najhNCtIZTorh//z6EQiFeeOEFpKSk9Jqpsd0p3IYPR1PxXUp+A9PpLChCCNEUTo+e1q5di7///e/49ddfsWPHDrz22mt44YUX8Morr2DYMNn1Gnql0kzAnntbk/G2BmwNxeqsZqU+I6mr0FQ9hQSX1uNUe0GI7uK8cBGPx4OdnR3s7OxgaGiIR48eYcaMGXj33XfVGR9RkmStCnWuU0EIIQDHRLFhwwYEBQXh3XffxciRI3Hjxg1s2bIF6enp+Omnn9QdI1GCaMRIPE54g13YiBBC1IXTo6fy8nL897//xTPPPCO1nc/n4/Dhw2oJjBBCSM/A6Y6iqalJJkmsXLkSAODh4dH1qPRFaSZQmgkv/j32ixBCdA2nRHH8+HGZbceOHetyMLrOoPkxBuZf0Mix206V1fSUWcPMDJkvQkjvpVKi2LJlC3x8fHD79m34+vqyX05OTr1q4SJ5tRTNffsBAKyK1P9Lte1UWQmaMksI6S4qjVG8/PLLmDhxIt577z2sX7+e3W5ubg4rKyu1B6dLmo0tYNRUByg3k1UlbafKSig7ZZYQQrpKpURhaWkJS0tL7N27V1Px9GxOo4C8sx2+nWPA4DVTERpyf5R5L6S/OyKsfTQZXY+mTO0F1VoQ0jOp9OgpPDwcgPgOwsLCgv2SvO7NQvq7w7WFJ/e9gsdlSKvKBiDuLGtamdWdoRFCSJeodEeRmipec6G2tlYjweiyCGsfLLr1OwDgduBMqfc+a3OH0bZpoNBO+i/o+nJxI8GsEmoiSAjpOThXZhNCCOkdVLqjMDc3B4/Hk9sEkMfjoaamRm2Bkc5Jpsy2tbhB/N9miwqD3eE2fIVrcxNCejeVEgU9cuo5wm34gBqmx+Y3MEB5KyUKQkiHVEoU4eHhSE1NZe8s2qM7ChWUZuJyaSaEAa+wm8xCRoi/OShb0NievCmzAGCSKv7v8pGnkVJhKDvNtqOiO1oNjxD9p9IYRdvB7JqaGpkvVaWkpMDNzQ0uLi5SdRntXb58GQYGBjhw4IDK5yCEENI1WlsKtaWlBUuXLsXx48chEAggFAoRGxsLT09Pmf1WrlyJ6OhoLUWqHgWPy/BZ7o/dUk8hWZdCQtPrU6gLrXNBSM/EadZTY2MjvvjiC0yfPh1xcXH48ssv0djYqNIxLl26BBcXFzg7O6NPnz6YNWsWkpOTZfb797//jbi4ONja2nIJVf2cRqn8kZD+7hhqMlCqnkJTJOtSSND6FISQruJ0RzFv3jyYm5vj9ddfBwDs3bsXc+fOxY8/ylYkd6S4uBhDhgxhXwsEAqSlpcnsc/DgQZw6dQqXL1/mEmqPEGHtgwhrH6l6Ck0RjRgpdffQ9s6CEEK44JQobt++jevXr7OvIyMj4efnp9IxOppi29Zbb72Fzz77DAYGimfkJCUlISkpCQBQVlamUhzqZlGRi4H5F1DmGKbcByQtQTjcqRBCSHfglCgCAgJw8eJFhIaGAgDS0tIwcqRqz8AFAgEKCwvZ10VFRRg8eLDUPleuXMGsWbMAiBdLOnr0KAwNDTF16lSp/RISEpCQkAAACA4OVvVy1KZSEAiLilxYFWUonyh6AHn1GG1RnQUhvZtKicLHxwc8Hg/Nzc3YtWsXhg4dCgAoKCiQGYTujFAoRE5ODvLy8uDg4IB9+/Zhz549Uvvk5eWx38fHx2Py5MkySaInKXMM00ibcU3qrB6D6iwIISolCnUuc2poaIiNGzciOjoaLS0tWLBgAby8vLB161YAQGJiotrOpQlt16S4XHVHrcdWVJug7kWEOqrHkKB25oQQlRJF2+VPHz16hJycHKnZTu2XR+1MTEwMYmJipLZ1lCB27typ0rEJIYSoB6cxim+++QZff/01ioqK4O/vj4sXLyIsLAynTp1Sd3w6SZUBbfZu5L4xu+12lbh7rFv/AI3Ep05t73C0UaWtbO0F1VsQwh2nRPH111/j8uXLCA0NxenTp5GdnY3Vq1erO7aeq/0MJckv+dLMTge0b9cX40zFDYVFd14OlrhZrL5W4+0L8NrSlWI8Qoj2cCq4MzY2hrGx+JdjU1MT3N3dcfv2bbUGpqvKHMNQYz1M7nsh/d0BQONFd221L8Bri4rxCCHK4HRHIRAIUFVVhalTp2L8+PEYMGCAzNRWIivC2qdbkwQgW4DXlrLFeFm10tNnabosIb0Lp0Rx8OBBAMCaNWsQGRmJ6upqTJgwQa2B6TqVC+96qHAbPoCn02dpuiwhvQ+nRNHY2IjNmzcjNTUVPB4P4eHhaG3t+toI+kJXC+/kaT99lqbLEtL7aK3Xkz4R2glx+f7TXlScCu9KM6VemlbWA3JmPbWdWaTumgpCCJFHa72eCCGE6AZOs54kvZ4kuPR6Ioo52ZhpOwRCCAGgxV5PvYFpdTHcUjdJbTMaYIxmYwstRdR7KSrMo2I8QhTTWq8nfVcpkK1SNq0uhpHpIOQyTd222l1n5BXjUREeIaQtzr2erl+/jrNnxWspjBo1isYo2ilzDJOZ8eSWugnjm1vRYCle7Q6AwkRxu+qqzDZ1tvVo9g+CUbtt/NJiGAEKE0WHbcmz0mQ2Mf36AQBGDrPBOI9BXYiWEKItnFt4bNu2DdOnTwcAzJkzBwkJCewsKNKxKc18uA+bqdRqd6aVWQCABivNPNaTV4zXWRFeZ23J5blX0QCgnBIFITqKU6LYvn070tLSYGYmHnBduXIlwsLCenWiENoJcbndFFdNa9+Erzumy3bWlrw9ka8XPj58U4MREUI0jdOsJ4ZhpJYnNTAwkLu0aa9j79vpLpIBbqPGmk73dbIxg5ONGbwcLNURHSGEcMLpjmL+/PkICQnBtGnTAACHDh3CK6+8otbA9JFkgFsyqE2znwghuoBToli2bBkiIiKQmpoKhmGwY8cOBAT0/LUTtE0ywC2eMivSdjiEEKIUlRMFwzAoKipCYGAgAgO7f6EaQggh3UvlMQoej4epU6dqIJTehS9qgml1CQbmX+h859JMmFZmwbQyCxb3L8Li/sXOP0MIIWrC6dFTaGgoLl++DKFQqO54eoVKQSBaKy+BL2rqkR1mNVGEd6+iodPZT9qqtVB2OVVFqLqb6DNOieL06dPYunUrHB0dYWZmBoZhwOPxkJnZvdNDdVWZYxgaWopgWl0CPFHfcSXTZbsyTZZrEZ4iI4fZAChXuA/VWhDSc3FKFMeOHVN3HKQTkiaBDbgPQPwXrDr+Em6PSxGeIoaZGYgGED1U8X6r61qBujqFSa593QghpHtwShSDBg2SWbho8eLF6o6NEEJID8Cp4G7evHm4efMmXn/9dbz22mu4desW5s6dq+7YegWLily4pW5SblCbEEK0gBYu0qLmvv1QYz0MFhW5AKDaoHbeWVjcF1d319iFaiI8QggBQAsXaVWzsQVuhy9FjfUwbYdCCCEd4nRHkZaWJrNwkYeHB7uwUa+e/WTvK7P+NSGE6DJOiSIlJUXdcegFod3/6koeN8q8d7nqjtrP5zX4f72inKTn8NeXq9ZEMKukWl0hdUmH61z8D1MgW4dB61wQonmcEkXbBYwId7fri3Gm4gbctB2IEuQV4QHqWw1P19e50MRU5a6gAkCiTpwSBem6kP7uuF1fjLSqbCzSdjCdkFeEB3S9EK8tZda5EPl6Sb2mdS4I6R6UKDTBaRSQd1bhLhHWPkiryu6mgLpGXhEe0LVCPEKI7uA060me+/fvq+tQRFWdJCVCCOkKtSUKWriIEEL0k9oePR05ckRdh+qVJBXabVUKAjsswpOZRXXfmP3WU8Vzew5WPEuqp8yKat8HilfXLHe7BPWGIkQ9ON1RrFy5UqltvZrTKKmXwv7DIew/XO6ulYJAmaI70+piWBVx6wJrFjJC6osQQrqC0x3F8ePH8dlnn0ltO3bsmMw2ohzJEqlttb+7IPIpqr1oW3dB9RaEcKfSHcWWLVvg4+OD7Oxs+Pr6sl9OTk7w8fFR+eQpKSlwc3ODi4sL1q9fL/P+7t272XM8++yzUv2lCAm34cPRlNfpfvcqGnAuV/F6GISQjql0R/Hyyy9j4sSJeO+996R+sZubm8PKykqlE7e0tGDp0qU4fvw4BAIBhEIhYmNj4en59Am7k5MTfvvtNwwYMADHjh1DQkIC0tLSVDpPT1fwuAyf5f4os93UVITxzXy4ayEmVbQtxFNX8Z2yOqu9kNRdUL0FIV2jUqKwtLSEpaUlpk+fDisrK5ibm2Pt2rXIyMjA3/72NwQEBCh9rEuXLsHFxQXOzs4AgFmzZiE5OVkqUTz77LPs96GhoSgqKlIlXO1rN04BALj6dBA6pH/HaSDHgAFf9BhTUjcpHNTWpraFeOosviOE9Cycxij+/ve/Y+bMmUhNTcUvv/yC5cuXIzExUaW/9ouLizFkyBD2tUAgUPj57du3Y+LEiXLfS0pKQlJSEgCgrKxM6Ri0LcLaBxHW8h/ZfXFzO1pb6mBaVgxAiRbkbRsRtu81df92V8KEZ7sHlPnl9eJvXAYALs+Jv9+XDL6oHqaVWR0ep8FK1flYhKue1lKEaJamW7ZwShQGBuLb/SNHjmDx4sWYMmUK1qxZo9IxGIaR2cbjyX/efPr0aWzfvh2pqaly309ISEBCQgIAIDg4WKU4upuw/3ClGgQ2G1ug2dgCDQ1dn8Fs5ttxN6n6zK4lEUKI/uM0PdbBwQGLFi3C/v37ERMTg6amJrS2qtbQTSAQoLCwkH1dVFSEwYMHy+yXmZmJhQsXIjk5GdbW1OiMEEK6G6dE8cMPPyA6Ohq//PIL+vfvj8rKSvzjH/9Q6RhCoRA5OTnIy8vDkydPsG/fPsTGxkrtU1BQgOnTp+P//u//MHy4/BoEneM0SrxmBSGE6AhOzzVMTExQX1+PvXv34sMPP0RzczP69++v2okNDbFx40ZER0ejpaUFCxYsgJeXF7Zu3QoASExMxMcff4yKigosWbKE/cyVK1e4hEx6uXsVDR3OfqIaC0IU45QolixZAj6fj1OnTuHDDz+Eubk54uLicPnyZZWOExMTg5iYGKltiYmJ7PfffPMNvvnmGy4h6hXT6mKV2nsQaSOH2QCQX0fRk9a0IKSn4rwUakZGBjsddsCAAXjy5IlaAyNilQLZfkWm1YpnQikaLO+ojYg+kvSAigYQPVT+PqvrWoG6ug77RfU01L+KaAOnRGFkZISWlhZ2llJZWRn4fLU1ou0dJOMUCtbXLnhchuUmAOxt2W0h/d2x6NbvGg5OMUcbM5ltD4wM8KSkDH1+Oiy13czfA/1CfJ9OqSWE6BxOieKNN97AtGnT8PDhQ3zwwQc4cOAA1q5dq+7YejV5xXgFj8U1Ij1xRTwzfw+ZbU9KxPH2CxEnRUU1FppCtRuEdB2nRDF79mwEBQXh5MmTYBgGhw4dgoeH7C8KIp/QTggAuHy/4zEdecV48lp9dJW8GgsutRX9QnzZhCDx4D/7OcdFCOk5OFdzubu7w929p3ci0gFKPIIihBBtojWzdVTbmVA0A4qQ3k3SskVTrTwoUeigtjOhOpsBRTqnaE0LZYXb8BV2siVEl1Gi0EFtFzrSlQWO5M2Uak8bM6PCbfhAuWrtZ9rLb2CA8lZKFERvUaIgvVpna1ooo6t3I4T0dJQoiAx1zYQihOgHShQ6RrIiXkh/9w7XsuhJnpSUKT1NltfcovB9xsMV8KO6CEK6GyUKHSIpwpMU3vX0RCGvCI+zhxXgAWBUTBRti/yo+I4QbihR9BTyWo+3q62QFOG1L7yzqMjFwPwLPW7mk7wiPEUUDWbz9iWrIyRCCAeUKHRcpSAQFhW5sCrK6HGJQlWKZkY9MBIPOA9qtw/1kCJE8yhR6LgyxzBYFelG51N9po5aDGUwBfLX1GiP1tgg6kQtXwnponAbPhxN5a/3rg33KhpwLlf++huEcEF3FFoktBMqbAzYk8ibMgvQtFlAPbUYyhL5enW6T0cr+RHCFd1R9GT2vrS+NiFE6+iOQk/IWy5VLkMTtZ7XuP6x0vuKBgRAZKPbA+6E9EaUKPSAvOVSexr+42IYAmpPFMr0kAJodhQhXUGJQkdJKrRZ7ZZL7agYT91rZjcqOUZhnLO5y+dSpcobeLoMKyGkayhR6IJ24xQhokrg3km5u+pK1baqVK3ybr8MKyGEO0oUWiZZFlWejmZERQyJQIShldz3NLFcqiJtZ0NpcgaUqlXevX0Z1nsVDTT7qRfRdN0MJQpC9MzIYTYAqI6it7hX0QCgnBIFIV3laGMG2FnKfS+rpLqbo9GscR6DqCq7F+mOO0dKFLqq7bhFu+aBPRX/cbFaBrWVOldDEwDAsPwCTcklpIuo4I50C9GAALSaOHTvSVubYPjoaveekxA9RHcUPZi8ge6e3PJD0cp4Ipuwbv3LvtV0P/gNJd12PkL0GSUK0ut5DpY/dqEKfRvnIKQtShR6SKYYr40XB4/BTPvwbo5IO5oqW1F4uAmtpv+bKtvniMbOZdIk0tix22qaNA1PJk7plnMRIkFjFPqgzcB2SH93DDUZKHe3gsdlOPqw5z66Uiczfw/0tdKvf9780mL0+e24tsMgPZCkbmZPWoFGjk93FHpGslyqPN1djKdN/UJ8YWOVCv7jYrSa9AUAiJ55FqKhkzRyvu549GSStEHj5yC6R1I3c6+iAcnXivFyyFC1n4MShY6RDHD35EHttjpax4ILVSu/RQMC2H/g/MfFMCw9pbFEQYi2SOpmNFlPQYlCXyizbkXhMaCPOeA0SvPxyJN3tltP13amlabrN9QxIN6ZB30NYWBsBJ9h1p3ueyG3QuPxkN5Dvx7iEkIIUTutJoqUlBS4ubnBxcUF69evl3mfYRi88cYbcHFxga+vLzIyMrQQJSGE9G5ae/TU0tKCpUuX4vjx4xAIBBAKhYiNjYWnpye7z7Fjx5CTk4OcnBykpaVh8eLFSEtL01bIPQqXrrNa18VHXmYqfr4+7RJwXzfam6hbmBKPp4j+0PSjRq3dUVy6dAkuLi5wdnZGnz59MGvWLCQnJ0vtk5ycjHnz5oHH4yE0NBRVVVUoLS3VUsSEENI7ae2Oori4GEOGDGFfCwQCmbsFefsUFxfD3t6+2+LURR3dbZj3Me/mSLTLLGQEgBHiFzsOttmmmwzMe9d/P6K8sGHWsDA20tjxtZYoGIaR2cbj8VTeBwCSkpKQlJQEAMjOzkZwcDCnmMrKyjBwoPxiNX0S/NenP5/ecs2sTcG6f80c/n3r/DVz0FuvOXgbt2vOz8/v8D2tJQqBQIDCwkL2dVFREQYPHqzyPgCQkJCAhISELscUHByMK1eudPk4uoSuuXega+4dNHXNWhujEAqFyMnJQV5eHp48eYJ9+/YhNjZWap/Y2Fjs2rULDMPg4sWLsLS0pMdOhBDSzbR2R2FoaIiNGzciOjoaLS0tWLBgAby8vLB161YAQGJiImJiYnD06FG4uLjA1NQUO3bs0Fa4hBDSa2m1MjsmJgYxMTFS2xITE9nveTweNm3a1G3xqOPxla6ha+4d6Jp7B01dM4+RN2JMCCGE/A+18CCEEKIQJYr/6aydiK4qLCxEZGQkPDw84OXlha+//hoAUFlZifHjx8PV1RXjx4/Ho0eP2M+sW7cOLi4ucHNzwy+//KKt0LukpaUFAQEBmDx5MgD9v14AqKqqwowZM+Du7g4PDw9cuHBBr6/7yy+/hJeXF7y9vfHSSy+hsbFRL693wYIFsLW1hbe3N7uNy3Wmp6fDx8cHLi4ueOONN+SWH3SIIYxIJGKcnZ2Z3NxcpqmpifH19WVu3ryp7bDUoqSkhElPT2cYhmFqamoYV1dX5ubNm8yKFSuYdevWMQzDMOvWrWPeffddhmEY5ubNm4yvry/T2NjI3L17l3F2dmZEIpHW4ufqX//6F/PSSy8xkyZNYhiG0fvrZRiGmTdvHrNt2zaGYRimqamJefTokd5ed1FREePo6Mg0NDQwDMMwM2fOZHbs2KGX1/vbb78x6enpjJeXF7uNy3UKhULm/PnzTGtrKzNhwgTm6NGjSsdAiYJhmPPnzzNRUVHs608//ZT59NNPtRiR5sTGxjK//vorM3z4cKakpIRhGHEyGT58OMMwstceFRXFnD9/XiuxclVYWMiMHTuWOXnyJJso9Pl6GYZhqqurGUdHR6a1tVVqu75ed1FRESMQCJiKigqmubmZmTRpEvPLL7/o7fXm5eVJJQpVr7OkpIRxc3Njt+/Zs4dJSEhQ+vz06AkdtwrRN/n5+bh69SpCQkLw4MEDtibF3t4eDx8+BKAfP4u33noLn3/+Ofj8p/+89fl6AeDu3bsYOHAg5s+fj4CAACxcuBD19fV6e90ODg5Yvnw5hg4dCnt7e1haWiIqKkpvr7c9Va+zuLgYAoFAZruyKFFA+VYhuqyurg5xcXH46quvYGFh0eF+uv6zOHz4MGxtbREUFKTU/rp+vRIikQgZGRlYvHgxrl69CjMzM4Vjbbp+3Y8ePUJycjLy8vJQUlKC+vp6fP/99x3ur+vXq6yOrrOr10+JAsq3CtFVzc3NiIuLw+zZszF9+nQAwKBBg9hOvKWlpbC1tQWg+z+Lc+fO4eeff4ajoyNmzZqFU6dOYc6cOXp7vRICgQACgQAhISEAgBkzZiAjI0Nvr/vEiRNwcnLCwIEDYWRkhOnTp+P8+fN6e73tqXqdAoEARUVFMtuVRYkCyrUT0VUMw+CVV16Bh4cHli1bxm6PjY3Fd999BwD47rvvMGXKFHb7vn370NTUhLy8POTk5GDECN3puLpu3ToUFRUhPz8f+/btw9ixY/H999/r7fVK2NnZYciQIbh9W7yu+MmTJ+Hp6am31z106FBcvHgRDQ0NYBgGJ0+ehIeHh95eb3uqXqe9vT3Mzc1x8eJFMAyDXbt2sZ9RCrehFf1z5MgRxtXVlXF2dmbWrl2r7XDU5uzZswwAxsfHh/Hz82P8/PyYI0eOMOXl5czYsWMZFxcXZuzYsUxFRQX7mbVr1zLOzs7M8OHDVZoZ0dOcPn2aHczuDdd79epVJigoiPHx8WGmTJnCVFZW6vV1f/jhh4ybmxvj5eXFzJkzh2lsbNTL6501axZjZ2fHGBoaMg4ODsw333zD6TovX77MeHl5Mc7OzszSpUtlJj4oQpXZhBBCFKJHT4QQQhSiREEIIUQhShSEEEIUokRBCCFEIUoUhBBCFKJEQQghRCFKFIQQQhSiREGIiqqqqrB582b29bPPPquR8xQVFWH//v0aOTYhqqBEQYiK2ieK8+fPa+Q8J0+eREZGhkaOTYgqKFEQoqJVq1YhNzcX/v7+WLFiBfr16wdA3Mbd3d0dCxcuhLe3N2bPno0TJ05g5MiRcHV1xaVLl9hjfP/99xgxYgT8/f2xaNEitLS0SJ0jNTUVy5Ytw4EDB+Dv74+8vLxuvUZCpKixJQkhvUL7RWTMzMzY7QYGBkxmZibT0tLCBAYGMvPnz2daW1uZQ4cOMVOmTGEYhmGysrKYyZMnM0+ePGEYhmEWL17MfPfddzLniY6OZm7cuKH5CyKkE4baTlSE6BMnJyf4+PgAALy8vDBu3DjweDz4+PggPz8fgPiRUnp6OoRCIQDg8ePHbJvotm7fvg03N7dui52QjlCiIESN+vbty37P5/PZ13w+HyKRCIC49ftf/vIXrFu3rsPjVFRUwNLSEkZGRpoNmBAl0BgFISoyNzdHbW0t58+PGzcOBw4cYJevrKysxL1796T2ycvL0+mFdYh+oURBiIqsra0xcuRIeHt7Y8WKFSp/3tPTE2vXrkVUVBR8fX0xfvx4drUyCXd3d5SXl8Pb21tjs6oIURatR0EIIUQhuqMghBCiECUKQgghClGiIIQQohAlCkIIIQpRoiCEEKIQJQpCCCEKUaIghBCiECUKQgghCv1/2wGnXkBmWC8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for value in data_x[\"Celltype\"].unique():\n",
    "    mask = data_x[\"Celltype\"] == value\n",
    "    time_cell, survival_prob_cell, conf_int = kaplan_meier_estimator(\n",
    "        data_y[\"Status\"][mask], data_y[\"Survival_in_days\"][mask], conf_type=\"log-log\"\n",
    "    )\n",
    "    plt.step(time_cell, survival_prob_cell, where=\"post\", label=f\"{value} (n = {mask.sum()})\")\n",
    "    plt.fill_between(time_cell, conf_int[0], conf_int[1], alpha=0.25, step=\"post\")\n",
    "\n",
    "plt.ylim(0, 1)\n",
    "plt.ylabel(r\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we observe a pronounced difference between two groups. Patients with *squamous* or *large* cells seem to have a better prognosis compared to patients with *small* or *adeno* cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Survival Models\n",
    "\n",
    "In the Kaplan-Meier approach used above, we estimated multiple survival curves by dividing the dataset into smaller sub-groups according to a variable. If we want to consider more than 1 or 2 variables, this approach quickly becomes infeasible, because subgroups will get very small. Instead, we can use a linear model, [Cox's proportional hazard's model](https://en.wikipedia.org/wiki/Proportional_hazards_model), to estimate the impact each variable has on survival.\n",
    "\n",
    "First however, we need to convert the categorical variables in the data set into numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Age_in_years</th>\n",
       "      <th>Celltype=large</th>\n",
       "      <th>Celltype=smallcell</th>\n",
       "      <th>Celltype=squamous</th>\n",
       "      <th>Karnofsky_score</th>\n",
       "      <th>Months_from_Diagnosis</th>\n",
       "      <th>Prior_therapy=yes</th>\n",
       "      <th>Treatment=test</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>69.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>64.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>38.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>63.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>65.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Age_in_years  Celltype=large  Celltype=smallcell  Celltype=squamous  \\\n",
       "0          69.0             0.0                 0.0                1.0   \n",
       "1          64.0             0.0                 0.0                1.0   \n",
       "2          38.0             0.0                 0.0                1.0   \n",
       "3          63.0             0.0                 0.0                1.0   \n",
       "4          65.0             0.0                 0.0                1.0   \n",
       "\n",
       "   Karnofsky_score  Months_from_Diagnosis  Prior_therapy=yes  Treatment=test  \n",
       "0             60.0                    7.0                0.0             0.0  \n",
       "1             70.0                    5.0                1.0             0.0  \n",
       "2             60.0                    3.0                0.0             0.0  \n",
       "3             60.0                    9.0                1.0             0.0  \n",
       "4             70.0                   11.0                1.0             0.0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.preprocessing import OneHotEncoder\n",
    "\n",
    "data_x_numeric = OneHotEncoder().fit_transform(data_x)\n",
    "data_x_numeric.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Survival models in **scikit-survival** follow the same rules as estimators in scikit-learn, i.e., they have a `fit` method, which expects a data matrix and a structured array of survival times and binary event indicators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CoxPHSurvivalAnalysis()"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import set_config\n",
    "\n",
    "from sksurv.linear_model import CoxPHSurvivalAnalysis\n",
    "\n",
    "set_config(display=\"text\")  # displays text representation of estimators\n",
    "\n",
    "estimator = CoxPHSurvivalAnalysis()\n",
    "estimator.fit(data_x_numeric, data_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a vector of coefficients, one for each variable, where each value corresponds to the [log hazard ratio](https://en.wikipedia.org/wiki/Hazard_ratio)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Age_in_years            -0.008549\n",
       "Celltype=large          -0.788672\n",
       "Celltype=smallcell      -0.331813\n",
       "Celltype=squamous       -1.188299\n",
       "Karnofsky_score         -0.032622\n",
       "Months_from_Diagnosis   -0.000092\n",
       "Prior_therapy=yes        0.072327\n",
       "Treatment=test           0.289936\n",
       "dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(estimator.coef_, index=data_x_numeric.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the fitted model, we can predict a patient-specific survival function, by passing an appropriate data matrix to the estimator's `predict_survival_function` method.\n",
    "\n",
    "First, let's create a set of four synthetic patients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Age_in_years</th>\n",
       "      <th>Celltype=large</th>\n",
       "      <th>Celltype=smallcell</th>\n",
       "      <th>Celltype=squamous</th>\n",
       "      <th>Karnofsky_score</th>\n",
       "      <th>Months_from_Diagnosis</th>\n",
       "      <th>Prior_therapy=yes</th>\n",
       "      <th>Treatment=test</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Age_in_years  Celltype=large  Celltype=smallcell  Celltype=squamous  \\\n",
       "1            65               0                   0                  1   \n",
       "2            65               0                   0                  1   \n",
       "3            65               0                   1                  0   \n",
       "4            65               0                   1                  0   \n",
       "\n",
       "   Karnofsky_score  Months_from_Diagnosis  Prior_therapy=yes  Treatment=test  \n",
       "1               60                      1                  0               1  \n",
       "2               60                      1                  0               0  \n",
       "3               60                      1                  0               0  \n",
       "4               60                      1                  0               1  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_new = pd.DataFrame.from_dict(\n",
    "    {\n",
    "        1: [65, 0, 0, 1, 60, 1, 0, 1],\n",
    "        2: [65, 0, 0, 1, 60, 1, 0, 0],\n",
    "        3: [65, 0, 1, 0, 60, 1, 0, 0],\n",
    "        4: [65, 0, 1, 0, 60, 1, 0, 1],\n",
    "    },\n",
    "    columns=data_x_numeric.columns,\n",
    "    orient=\"index\",\n",
    ")\n",
    "x_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to `kaplan_meier_estimator`, the `predict_survival_function` method returns a sequence of step functions, which we can plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x121a78b90>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAENCAYAAAARyyJwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAAsTAAALEwEAmpwYAABR20lEQVR4nO3dd1hU19bA4d8MHaSIFIFBEUFEBFFBscRYYkk01sSYGBPTjBpTruk3N12jKfcmfjHRmGJM0xiNGnuLHQvYiA1RUSkqXen1fH+MokiRGWao632eeWDO2WefdcaENafstVWKoigIIYQQlVDXdQBCCCHqN0kUQgghqiSJQgghRJUkUQghhKiSJAohhBBVkkQhhBCiSqZ1HYChOTk54eXlVddhCCFEg3L+/HlSUlIqXNfoEoWXlxeRkZF1HYYQQjQoISEhla6TS09CCCGqJIlCCCFElSRRCCGEqFKd3qN48sknWbNmDS4uLhw7dqzcekVRePHFF1m3bh3W1tb8+OOPdOnSpQ4iFULUB4WFhcTHx5OXl1fXoTRYlpaWaDQazMzMqr1NnSaKiRMnMm3aNB577LEK169fv56YmBhiYmLYv38/U6ZMYf/+/bUcpRCivoiPj8fW1hYvLy9UKlVdh9PgKIpCamoq8fHxtGnTptrb1emlpz59+uDo6Fjp+lWrVvHYY4+hUqkICwsjIyODS5cu1WKEQoj6JC8vjxYtWkiS0JNKpaJFixY6n5HV63sUCQkJeHp6lr7XaDQkJCTUYURCiLomSaJm9Pn86nWiqGiqjIoOcsGCBYSEhBASEkJycrJe+8rLSifu+wmk7f5er+2FEE3DzJkzCQgIICgoiODgYKNfDu/bt69OY8P++OMPAgICUKvVBhtTVq8H3Gk0GuLi4krfx8fH4+7uXq7dpEmTmDRpElD1oJGqpGVlsPHadnruXE9zK1NUnSeAul7nUSFELdu7dy9r1qzh0KFDWFhYkJKSQkFBQV2HVUbHjh35888/efbZZw3WZ73+Szh8+HB++uknFEVh37592Nvb4+bmZpR9qRVwW9mMXxNdUK1+gbyLB42yHyFEw3Xp0iWcnJywsLAAtCWDbnx5/eCDDwgNDaVjx45MmjSp9IpI3759+de//kWfPn3w9/cnIiKC0aNH4+vry3/+8x9AWz6jffv2PP744wQFBfHAAw+Qk5NTbv+bNm2iR48edOnShQcffJCsrKxybfz9/fHz8zPocddponj44Yfp0aMH0dHRaDQavv/+e+bPn8/8+fMBuO+++/D29sbHx4dnnnmGr7/+2mixONm0wPsKPLBDIVOlIi1iqdH2JYRomAYNGkRcXBzt2rVj6tSp7Nixo3TdtGnTiIiI4NixY+Tm5rJmzZrSdebm5uzcuZPJkyczYsQIvvrqK44dO8aPP/5IamoqANHR0UyaNImoqCjs7OzK/b1LSUlhxowZbNmyhUOHDhESEsL//ve/WjnuOr30tHjx4irXq1Qqvvrqq1qJxdTOjjwXeyyTrjKzuSPj0tMpf5FLCFFfvL/6OCcSrxm0zw7udrx7f0Cl65s1a8bBgwfZtWsX27Zt46GHHmL27NlMnDiRbdu28cknn5CTk0NaWhoBAQHcf//9gPbqCEBgYCABAQGlV0a8vb2Ji4vDwcEBT09PevXqBcCjjz7K//3f//HKK6+U7nvfvn2cOHGitE1BQQE9evQw6PFXpl7fo6htzUeOInfBj8RmWKEyWwXZs8HGqa7DEkLUIyYmJvTt25e+ffsSGBjIokWLGDduHFOnTiUyMhJPT0/ee++9Mo+g3rhUpVarS3+/8b6oqAgo/6DO7e8VRWHgwIF3/IJtDJIobuFyV38uLPiRV5eXUDBaDWtfhrGL6josIUQFqvrmbyzR0dGo1Wp8fX0BOHLkCK1bty5NCk5OTmRlZbFs2TIeeOABnfq+ePEie/fupUePHixevJjevXuXWR8WFsZzzz3HmTNn8PHxIScnh/j4eNq1a2eYg6tCvb6ZXdssAwJQbKywyYeSDFM4sbKuQxJC1CNZWVk8/vjjdOjQgaCgIE6cOMF7772Hg4MDzzzzDIGBgYwcOZLQ0FCd+/b392fRokUEBQWRlpbGlClTyqx3dnbmxx9/5OGHHyYoKIiwsDBOnTpVrp8VK1ag0WjYu3cvQ4cOZfDgwXof7w0qpaLBCg1YSEhIjZ4dPrZpCSYvvM/KESa8aRUHE9eCV+87byiEMLqTJ0/i7+9f12EY3Pnz5xk2bFiFNe+MoaLPsaq/nXJGcRsvOy8ALpRotAuWPVV3wQghRD0gieI2KrQ3kJrlWJOh2EDWZThS+zePhBBNh5eXV62dTehDEsVtzDy0D8X6Xs3lkYK3tAtPrKrDiIQQom5JoriNuacnueZwTXWBE6qWFGAOp9fXdVhCCFFnJFFUQK0yAUDjlE9ksY92YUZcFVsIIUTjJYmiAmYm2uElHb0KWVfSXbtQHpUVQjRRkigqoLr+sVg1P8IW9fUh8qfW1mFEQoj6or6XGX/11Vdp3749QUFBjBo1ioyMjBrHIImiAuqiYoZFKKSkxZFSbKtdeHEvXNhbt4EJIerUrWXGo6Ki2LJlS5nJ1eqDgQMHcuzYMaKiomjXrh2zZs2qcZ+SKCpg4d8eANOkDPr6ufBcwQvaFZE/1GFUQoi61hDKjA8aNAhTU+3l87CwMOLj42t83JIoKtBi4kQALuVcprWTCZtKrk+GdLXmH7gQouFqaGXGf/jhB+69994aH7cUBazI9W8CfvEK3e6z4/tdpsQ1C8LzYjgU5oKZVR0HKIRg/Rtw+R/D9tkyEO6dXenqhlRmfObMmZiamjJ+/PiafSZIoqiQVadOgDZRnMpZC7Qlq/D6yp9Hw5MyrkKIpqohlBlftGgRa9asYevWreX60YckigqYuroC0OaKwvqrZ4C2jL36Av9YPg0XwyH9AjRvXbdBCtHUVfHN31gaQpnxDRs28PHHH7Njxw6sra1rcLQ3yT2KCqhMTbHq1Am1pRV7E/cyrZ8PmVgT3v56SY+4A3UboBCiTjSEMuPTpk0jMzOTgQMHEhwczOTJk/U+3hukzHglLj75FJl7wxn3himr7t9E/08OMdIzmy+Sn4Euj8Pw/zNAtEIIXUiZccOQMuMGUpSWhloB6zyFJacXApBk3kq78tAiyEmrw+iEEKL2SKKohMPo0QCYlqgoKCkgpHVzMnIKoe0AbYMvu9RhdEKIxkTKjDdUau1HE5BpB0BRicKJS9fY0ulz7frcdNjwbzi1rq4iFEKIWlGjRJGdnU1xcbGhYqlXzNy1zzn335PF+tj1vHiP9imHp387Bg/9qm207ytY8jDklx8dKYQQjYVOiaKkpITffvuNoUOH4uLiQvv27XFzcyMgIIBXX32VmJgYY8VZ62z798e8dWsUwNzEnH5+Lng72QBwzO4ueCcdel4v7bHv68o7EkKIBk6nRNGvXz/Onj3LrFmzuHz5MnFxcSQlJbFr1y7CwsJ44403+OWXX4wVa61TWVrS/pKajPwMFEXhraHapwSSMvO0l6b6XB81uW2m9lKUEEI0Qjolii1btvD2229jb2+PWn1zU0dHR8aMGcPy5ct56KGHDB5kXSlKTsbymnYgTezVWFztLAHYcOyytoGlPbTQXpLi6551EaIQopbV9zLjb7/9dmlsgwYNIjExscYx6JQozMzMABg1alS5dfv27SvTpjFweGgsAAOOlJBdmE2Au/bGdmxK9s1GT2/R/sxMhMTDtR2iEKIWNYQy46+++ipRUVEcOXKEYcOG8cEHH9S4T50SxdKlS3njjTfIzMzk5MmTZW5kT5o0qcbB1Df292sLebW9pBCTEVNaMyXifDoHL1wfR2HlAA9ox1mw+/M6iFIIUVsaQplxOzu70t+zs7MNUutJp0TRq1cvOnToQHp6OtOnT8fX15cuXbowbNgwrKwaX0VVC+82qJwcAYjL1M6Z/ckDQQBM/CHiZsOOo6G5V22HJ4SoZQ2lzPhbb72Fp6cnv/76q0HOKHQqCujh4cFjjz1G27ZtS0vdpqWlERsbS/v27WscTH1kojJBjRoFbVZ+sKuGd1cdJzO/iE3HLxPgYY+HgxWYWMCZrXUcrRBNx8cHPuZUWvlaRzXR3rE9r3d7vdL1DaXM+MyZM5k5cyazZs1i7ty5vP/++zX6XHQ6o7hxKnUjUNDeyO7atSs2NjZl2jQWSkkJ7snFrDq7CtCW/v339aefJv18kH6fbdc2zMuAAhlPIURjd6PM+Pvvv8/cuXNZvnw5eXl5TJ06lWXLlvHPP//wzDPPGK3M+JEjRzhy5AgnTpzg+++/rzLWRx55hOXLl9foeEHHM4p+/foxZswYRowYQatWrUqXFxQUsHv3bhYtWkS/fv2YeH2GuMagODWVNmaQlJNESm4KTlZOPBzqiX9LW15ccoSEjFxOXrqGf/AjsPerug5XiCajqm/+xtIQyozHxMSUxvfXX38Z5GqPTmcUGzZswMTEhIcffhh3d3c6dOhAmzZt8PX1ZfHixfzrX//SKUls2LABPz8/fHx8mD27fG35q1evcv/999OpUycCAgJYuHChLuEahN3QoXD9G8Cm85sAMDVRE+LlyIyRHQFYFH5e27i4AArzKupGCNEINIQy42+88QYdO3YkKCiITZs2MWfOHL2Pt5Sip4KCAiUxMVFJT0/Xa/uioiLF29tbOXv2rJKfn68EBQUpx48fL9Nm5syZymuvvaYoiqIkJSUpzZs3V/Lz86vst2vXrnrFU5nEd95VTnbqpHT8saPy8/Gfyx5DcYnS+vU1SuvX1yjFG99RlHftFOXo7wbdvxDiphMnTtR1CEYRGxurBAQE1Nr+Kvocq/rbWeUZxYkTJ3j00UdL3w8YMIDjx48D2vESERERzJ07lwMHdJ/I58CBA/j4+ODt7Y25uTnjxo1j1apVZdqoVCoyMzNRFIWsrCwcHR0xNa3dSfmUvFyUvHxc0xXCE8PLrDNRq/B30z6K9n7C9Wqy57bXanxCCGFsVSaKAQMGMGPGjNL38fHxBAQEABAeHs6ECRO4ePEiEydOZMWKFTrtOCEhocxAFY1GQ0JCQpk206ZN4+TJk7i7uxMYGMicOXPKjAivDdbdugPQLBfyi/PLrV8woSsAi6JNUdSmkHmpVuMTQjR8DbrM+KZNm3jrrbdK3986kOOnn35i8uTJLFiwgO3bt/Pxxx/rtGOlgqejbr/Lv3HjRoKDg0lMTOTIkSNMmzaNa9euldtuwYIFhISEEBISQnJysk5x3ImpUwsA2jf3Q0X5gSuejtaM7uIBQLGihrN/w+X6+w8uhBC6qjJRBAYG8uuvv5a+9/HxYdmyZSQlJbFy5UpGjBgBgIuLC/n55b9tV0Wj0RAXF1f6Pj4+vnSE4w0LFy5k9OjRqFQqfHx8aNOmTYU3byZNmkRkZCSRkZE4OzvrFEd1eZ7LrHTd8/21TxisaPG0dsH8XpBfeXshhGhIdLqO8/nnn/PNN9/g4eFBly5d6NlTWwivsLCwwqHkVQkNDSUmJobY2FgKCgpYsmRJ6aCUG1q1asXWrdpBbFeuXCE6Ohpvb2+d9lNTlkHakdgWOUVEXomksLiwXJs2Tjb4ujRjl/1w8LpLu3B97T+6J4QQxqBTomjZsiWbN28mPz+fdetuzuy2bds2+vXrp9OOTU1NmTt3LoMHD8bf35+xY8cSEBDA/PnzmT9/PqCtghgeHk5gYCADBgzg448/xsnJSaf91JRp8+ZgagooFCvF7EzYWWE7BbhwrQTG/qRdcORXSDlTa3EKIYSxqJSKbhY0YCEhITqV5K2Ok+39oUcXxvaNoqVNSzY/sLn8fmdsISUrn/Ozh8KW92H3/8C9C9i5g70GhswGAxTnEqIpO3nyJP7+/nUaw8yZM/ntt98wMTFBrVbzzTff0L17d6Ptr2/fvnz22WeEhITotN1nn33Gq6++SnJycrkv2BV9jlX97azdZ00bMLOEFAAuZ1+mRClBrSp7MnZ3O2eWH4rXvrlrOsRHaCczio+EU2ug+7PgWLuXzYQQhnVrmXELCwtSUlIoKCio67DKiYuLY/PmzWUqaNRE7T5r2kA1u/tuCi9e5EmH+wCISo4q18bFzgJzk+sfp4UtTFwDU/bAgHe0y/Z+DSfXyE1uIRqwhlBmHOBf//oXn3zyiUFKjIMkimqx6qIdTNf3YjMA8oorLtNRUFxCUXFJ2YV+94LKBCK+hd/Hw/5vjBqrEMJ4GkKZ8b/++gsPDw86depksOPW6dKTra1thRlKURRUKlWFYxwag+bjHiL5889LB/vlF5V/FLigSJsgIi+kE+bd4uYKa0d46R/tZaj5vaFIakEJYQiXP/qI/JOGLTNu4d+elv/+d6Xr63uZ8ZycHGbOnMmmTZsM96GgY6LIzGzal03MS7SJYuvFrdzteXeZdcOC3Ph+dyyrjiSWTRQA9h7aFwrs/BR6Twdz61qKWghhSDfKjPft25fAwEAWLVrEuHHjmDp1KpGRkXh6evLee+8Zrcz44sWLK43t7NmzxMbGlp5NxMfH06VLFw4cOEDLli31Pma9b2anp6cTExNT5sPo06eP3oHUZ6rr9aXMN4bDg2BtVv6PfEcPewBM1VVcE/ToCgkH4dgy6PKYUWIVoqmo6pu/sdT3MuOBgYEkJSWVvvfy8iIyMrLGwwr0ShTfffcdc+bMIT4+nuDgYPbt20ePHj34+++/axRMfaW2sUFlbY2pkxNWpldIyEwo18bMRI2dpSmRF9JLL8WVc+8n8N0A2PgWdBwD5ja1EL0QwlCysrJ4/vnnycjIwNTUFB8fHxYsWFCmzLiXl1eNyow/++yz+Pr6Vllm/EYljBkzZpSbj8IY9BpHERgYSEREBGFhYRw5coRTp07x7rvv8vvvvxsjRp0YYxwFwPnxj6IyM2N4/0OUKCVsGrMJt2ZuZdr4/HsdRSUKW1++m7bOzSru6KeRcG4bhDwJwz43eJxCNGb1YRyFMZw/f55hw4bVWmFAXcdR6PXUk6WlJZaWlgDk5+fTvn17oqOj9emqwZnSSZvlF0eXv044c5R2IqMhX+wkr7C44g5GXX/qKfIHOLGq4jZCCFGP6JUoNBoNGRkZjBw5koEDBzJixIhyBf0aq3vb3AvAtovbyq8LdMPVzoLCYoX1xyopN27rCvf/n/b3da9C4xoYL4TQQ4MuM16ZFStW4ODgwHvvvceHH37IU089xcqVKw0cWj2jKOTs24fbVTVdXLpw/tp5cgrLDoixszTj6/HaMRdz/66izlPXx0GlhqwrkCr1oIQQ9ZteieLzzz8nPl5bruLuu+9m+PDhmJubGzSw+sbq+uNmuUej6OSi/f3jiI/56+xfZdp1be1IBzc7ziZnE342pfIOx3yv/fl1mMyzLYQOGll5ulqnz+enV6K4du0agwcP5q677uKrr77iypUr+nTToDiMfRCAyx9+SB+PPjhaOrL67Go+i/isXNuX7tE+Ovf+Xyd45Nt9vPLHUUpKbvvH8RkAtu5QUiRlPYSoJktLS1JTUyVZ6ElRFFJTU0vvMVdXjarHRkVF8fvvv7N8+XI0Gg1btmzRtyuDMdZTT0pxMac6BgLgf/IEAO/vfZ9lp5exccxG3JvdvEeTW1DM84sPczW3gMvX8ohLyyXirXtwtrUo2+mBb2HdK/DKGWhmnAmXhGhMCgsLiY+PLzN+S+jG0tISjUaDmZlZmeVGqx7r4uJCy5YtadGiRZlBHo2RysSE5uPHk/7LL+SfPYtF27b4OmjPHI6lHCuTKKzMTfjucW1J4IV7Ynl/9Ql2xSQzuoumbKdqE+3POZ3gqU3QsmOtHIsQDZWZmRlt2rSp6zCaHL0uPc2bN4++ffsyYMAAUlJS+Pbbb4mKKl9RtbG5cZ+i6HpSDGmpTQaHkw5Xus2A9q4AfLTuFAP+u53t0bckVP/hEPo0FGZDxkUjRS2EEDWjV6K4cOECX3zxBcePH+f999+nQ4cOho6rXjJzK1srpV1z7YjITecrL8Dl0dyKJ3p50aNtC84mZ3PwQvrNlTZO0HmCUWIVQghD0evS0+zZsw0dR4NSfO3mzWcLEwuScpO4lHWp3EhtABO1infvDwBgbVRi5Z2e2Qzt7zN4rEIIUVM6nVHcKFJla2uLnZ1d6evG+8ZObas9xqxdN+fNHt5WWz74/X3v33H7EgW+/PsMuQW3jNp2bq/9eX43VFC+XAgh6ppOiWL37t2Attz4tWvXSl833jd2ln7aS01qq5vVY9/s9iaOlo7sSdhzx0f2gj0dADh5+ZbPyswSnNpBymkI/9LgMQshRE3pPeAuIaF8BdWmQH3bmZOZiRne9tq5sLMKK56W8IbJd7cF4LONt9XFmrBC+3P7bCi5bYY8IYSoY3oPuBs0aFCTGnB3g1JURPrPP6MUFpYuG9BqAADf/vNtldsODnDF0cac8LOp7DydfHOFvQbMbKCkEDIuGCVuIYTQl16J4t133+X48eN89dVXJCYmcvfdd3PPPfcYOrZ6ycReO0FR/rnY0mWDvQYDEJ8ZX+W2KpWKp3prnwFfGhlXduWNkuPfDYB98wwUrRBC1JxeieKGpjTg7gbXN94AIO2nRaXLnK2daWvftlrbP9fPB28nG9ZEXeLnvedvrmjbD7o9C8WFcHJNpdsLIURtkwF3OrIK1g66y48+rXcfH4zQjsBeEnHLWUUzF7jvEzCzggu7axSjEEIYks7jKBRFITIyki+++ILg4GAjhFS/mbm6YtOrF9l79lCSn4/awuLOG92mt692/trjide4mJpDqxa3zMHduqdMaCSEqFd0PqNQqVQcPny4SSaJG0xdXAAovHgRpaAAAAWFLReqXxTxxmx4c7bGlH2s1qEVmDTuku1CiIZFr0tPPXr0ICIiwtCxNBjWXbWTE527fzixYx8CIL84H4XqF+Lt7NkcgOWH4knPufkEFYoCRXlwteob40IIUVv0ShTbtm0jLCyMtm3bEhQURGBgIEFBQYaOrd6yHTKElu+/j1VIV/JPnQJgiNcQzNRmd9jypg7udnw4Qlvao+TWMwp7T+3P2J0yoZEQol7Qq9bT+vXrDR1Hg2LSrBnNHxpL7uHD5EYeJD8mpkb9lUkUGm1FWlZOgWN/wqPLatS3EELUlF6JYtGiRRUuf+edd2oUTEPTrG9frq5cSUF8POhR6kqtVgHQ++NtrJjakwB3e3ALhnG/wbZZkJ1cdQdCCFEL9Lr0ZGNjU/oyMTFh/fr1nD9/Xud+NmzYgJ+fHz4+PpVWpN2+fTvBwcEEBARw99136xOu0Zh5aiciyt67F4DCksKqmpdzX0c3nujlRUFRCYkZ1y8zqdXQfijYuVe9sRBC1BK9zihefvnlMu9feeUVhg8frlMfxcXFPPfcc2zevBmNRkNoaCjDhw8vM7dFRkYGU6dOZcOGDbRq1areDeqzCtDeY8g7dpz8e7T3aE6lnaK9Y/tqbd/cxpwxXTQs3HO+4gaXjkBxEZjUaCJCIYSokRqNzL4hJyeHc+fO6bTNgQMH8PHxwdvbG3Nzc8aNG8eqVWXHD/z222+MHj2aVq1aAdqR4PVRUUpK6Wx3STn6JbOMnIKyC0yu3xjPrl/JUQjR9OiVKG485RQUFERAQAB+fn68+OKLOvWRkJCAp6dn6XuNRlOuIu3p06dJT0+nb9++dO3alZ9++kmfcI3KdsgQCi9exD38LAAbYjfotr2l9mxh++nb7ke009aP4v86Q2LlU60KIYSx6XVNY82am7WITE1NcXV1xdRUt64qmrtBpVKVeV9UVMTBgwfZunUrubm59OjRg7CwMNq1a1em3YIFC1iwYAEAycm1ewO42V13kblhA1Zrd8M9YG1mfeeNbtG6hQ1qFWw8dpnkzHycba+P9G4/DK6cgP3z4PRGcPIDc936FkIIQ9DrjOLAgQM4OjrSunVrFi5cyNixYzl06JBOfWg0GuLibtY6io+Px93dvVybIUOGYGNjg5OTE3369OHo0aPl+po0aRKRkZFERkbi7OyszyHpzWHMaKxDQsiNiMDRzEGvPoI0DhSVKGWLBFo7Qthk7e/bZ8mkRkKIOqNXovjwww+xtbVl9+7dbNy4kccff5wpU6bo1EdoaCgxMTHExsZSUFDAkiVLyt0QHzFiBLt27aKoqIicnBz279+Pv7+/PiEblep6vSfT/CLCE8PvONPd7ZY+2wOAxRFxZbdt7gVT92tLesTth6O/Q0GOocIWQohq0StRmJiYALB27VqmTJnCiBEjKCgouMNWZZmamjJ37lwGDx6Mv78/Y8eOJSAggPnz5zN//nwA/P39GTJkCEFBQXTr1o2nn36ajh076hOyUdncpZ1LvPPRLOIy44i9GnuHLcoyN9X+MyRn5hOfnlt2pUt7sHWDs1thxSSIXmeQmIUQorpUiq5ff4Fhw4bh4eHBli1bOHjwIFZWVnTr1q3Cy0K1LSQkhMjIyFrdZ150NLEjRpLVK5An+5zkUf9Heb3b6zr18eeheKYvPcpjPVrzUKindvDdDQXZcPkf+GEw9JgGg2ca+AiEEE1dVX879TqjWLp0KYMHD2bDhg04ODiQlpbGp59+WqMgGzJLPz/MfdpiY24DwC8nfyElN0WnPlq3sMbSTM1Pey/wzY7bHjU2twGX6+NLUs8aImQhhKg2vRKFtbU1o0ePxtfXFwA3NzcGDRpk0MAaIjtze/7d/d8A7E3cq9O2XVs7curDe/F2smFnTAVPblnaaZNF0glY+wqk1Ky+lBBCVJdBBtyJm/wdtTfb15zTbzrT9JwCMnIqKQXSKgwKcyDiW/hHigUKIWqHJAoDC3YJpmOLjoQnhpOam6rz9mNDPbEwreSfZdjn8OoZ7e8xG2sQpRBCVJ9OiWLChAkAzJkzxyjBNGgKZP39NwCettoR5wevHDTe/oqLjNe3EELcQqdEcfDgQS5cuMAPP/xAeno6aWlpZV5NWUlODkqh9pLRM0HPAPDdP98xbes0ZuybodPYivyiEq5WdvkJwO8+uPIPHP61RjELIUR16JQoJk+ezJAhQzh16hRdu3Yt8woJCTFWjA2C/bBhqMy0hfw8mnkQ5hYGaKvJ/h79O9cKrlWrHztLbR9RCRmVN+owQvvzmNynEEIYn06J4oUXXuDkyZM8+eSTnDt3jtjY2NKXrtVjGzNrM2u+HfQtS+9fyhMdnwDgwrUL1dq2extHAD7dGM0v+yrZptM48AiBs39r59gWQggj0utm9rx58zh69Chz585l7ty5REVFGTquRsOjmQcAZzOqN/6hrXMzevs4EZuSze8RcXfeoCi/JuEJIcQd6ZUo/u///o/x48eTlJREUlIS48eP58svpWhdRXyba8eazDs6j4fWPMS+S/uqbN/cxpxfnu5OqJdj1R23H6r9+fMo+H0CFOZW3V4IIfSkV5nx7777jv3792Njox2J/Prrr9OjRw+ef/55gwbX0CiFhSiKUqZcuqu1K2N8x5Cam8r2+O1EXI4ovX9RI953w7k+kHkZLoZD+lvaulBCCGFgep1RKIpSWhgQtEUC9SgZ1aiU5GnnvM49fKTMclO1Ke/1fI8vB3yJWlX9j7tEUfgn4Sr5RcUVN/DoCo+vhr5vat9HLdEnbCGEuCO9ziieeOIJunfvzqhRowBYuXIlTz31lEEDa2hsBwwg/eefKcnOMkh/N55+ysgpxNXOpPKGbp20P8/vNsh+hRDidnqdUUyfPp2FCxfi6OhI8+bNWbhwIS+99JKBQ2tY1Ncvw11bX/lUqCVKCQuiFpBbdOf7CWHeLQB4YH44Z5IyK2/Yoi1494P4CIg7oFvQQghRDXqdUQB06dKFLl26GDKWBs2yg7bGk9rSotI2QU5BRKVEkZ6XjlUzqyr769femWFBbqyJusSZpCx8XGwrb+zWCc5tgz1zYJwMwhNCGJbUejIQlYkJJs2bV9nmgXYPALAk+s73E9zsrZjW3weAl5ceZdWRhMobD3wfXAPhWiKcXA151RvcJ4QQ1SGJohaFtgwF4GTqyWq193Fuxkv3+JJXVMKJxDv88bdygMRD8PujsP+bGkYqhBA36ZUo5s6dS3p6uqFjafCU4mLSf1tMcWbF9xQ0tho6OXdi36V97Irfdcf+TE3UvHRPO8xMVPweGcd9c3Zx6nIlCePhxTAlHFQmUI17IEIIUV16JYrLly8TGhrK2LFj2bBhQ5N/NPYGy3btACiMq3xEdf9W/QHd5qt4vr8vwZ4OnLh0jWMJlSQKC1twDYBbxnAIIYQh6JUoZsyYQUxMDE899RQ//vgjvr6+/Pvf/+bs2aY9TafjExMBuPT++yjFFY9/eLLjk2iaaVgXu44z6Weq1e9z/Xz4cERHAHZVNPudEEIYkd73KFQqFS1btqRly5aYmpqSnp7OAw88wGuvvWbI+BoUCz/tyOi8o1EUpVQ+Z3awSzAAy2OWV7tvVztLbd+FlQzAu6GkCHb9F/INM55DCCH0rvXUtWtXXnvtNXr16sU///zDvHnzOHjwIMuXV/+PX2NjrvGg5QfvA5C9e0+l7WbdNQsrUyv+vvh3tS/bmZuq8XO15VjCNV5fFsXpK5WMrWjVU/szp/JEJYQQutBrHEVKSgp//vknrVu3LrNcrVazZo1+c0U3FpZ+fgBcW7cOhzGjK22XW5RLblEuF65dwMveq1p99/RpwabjV/g9Mg5Xe0umD6xgbEWXCdraT2umg7m1dpmtGwz5GNTykJsQQnd6/eXIz88vlyRef/11APz9/WseVQNm1akTloGBZO/ZQ0lOTqXtPur9EQAJWVWMj7jNu/cHsOcN7c3wQxcqeerMLVj7yrwEqWchLgIOLIBsubchhNCPXoli8+bN5ZatX7++xsE0FmpL7f2E7PDwStvYW9gDMHnLZHIKK08olUnKzKt4hWsHeHYHTN2rfd19/Z5R1BKIXAiJR3TelxCiadMpUcybN4/AwECio6MJCgoqfbVp04agoCBjxdjguP7nPwBk79vP1dWryd5fvgZTT/ee9NX0BahW7adb3ePvyukrWZxNrsYNa9uW2p+b34E1L8Gq53TalxBC6HSP4pFHHuHee+/lzTffZPbs2aXLbW1tcXS8w0Q7TYiJgz2o1aT/8gvpvwAqFe0iIjBpZlPaxlRtSi+PXmyP365z/y522npSH609yfcTQ6tu3H4ovHoWigthzb8g9QyU3PrklEruXQghqqRTorC3t8fe3p7FixcbK55GwczVFd+dOyjJyiJj5UpS53+DUlgA2Nxx2+qYMaIjO08nE59ezTMRGyftT3MbSI2BD25J6jbO8OJR7TohhKiATl8le/fuDWjPIOzs7EpfN96Lm0ydnDD38sL0eqHA7N2Vzxfx0raXSMxKrHbfarUKc1M10VcyKSwuqX5Qd02Hfm/dfPndp73JnZtR/T6EEE2OTmcUu6//scuspJaRKM+mTx+YNZuiK1coiE/AzMO9dKrUbm7dCG0ZSsTlCKLTonFv5l7tfvv4OnMuOZv07ALsrMywNKticqMbXAO0rxsOLoLoddonpOw9dD00IUQTIRenjezGhEZJn/2Xs/fcQ/ovN+eL8Lb35pWQV/Tq18NBO59Ft4+2EvzBJlKy8nXvxNZN+zPhkF4xCCGaBp3OKGxtbVGpVBWOJlapVFy7JvMg3M7MxQXPbxdQlJrKpTfepCgt1SD9jumqwcJMzdG4qyw/FE96dgFOzSqfNKlC7p21P3NSIOWWulM2Ttqy5UIIgY5nFJmZmVy7do3MzMxyL32SxIYNG/Dz88PHx6fMU1S3i4iIwMTEhGXLlum8j/qg2V134TByZJVPF2UW6nY5z9HGnMd6eNGvvTMAX2/XoyCjiXZebnZ8DHO73vIK0b0vIUSjpdMZRe/evdm9e3fpmcXtdEkWxcXFPPfcc2zevBmNRkNoaCjDhw+nQ4cO5dq9/vrrDB48WJdQGwxbM20Zjp3xOxnedrjO23dvo51b+3xqtu47t3KAx1dD5pWby06s1N63EEKI6+rsZvaBAwfw8fHB29sbgHHjxrFq1apyieLLL79kzJgxRERE1Hifda6khNR583F6+unSexeedp6oUGGq1m/6cmdbC+7ydWJXTAoHYtPo1kbH8Sxt+pR9nxINp9bATyNvLuv2jHY8hhCiSaqzm9kJCQl4enqWvtdoNCQkJJRrs2LFCiZPnlzb4RmFdbduABRnZJRZ7mrjytpza1l2Wr9La8GeDgB8t+tcTcLT8u4HrXpAQbb2dSEcjv1Z836FEA2WXokiLy+P//3vf4wePZoxY8bw+eefk5dXSe2hSlR2Q/xWL730Eh9//DEmJlU/+rlgwQJCQkIICQkhObn+Fr+zHzECgIxVq8osf7+ntjR5fGa8Xv2+PMgPfzc7MnIKORKXodvYitt59YInN8DTm7UvO3e4dASKCvTvUwjRoOmVKB577DGOHz/O888/z7Rp0zh58iQTJkzQqQ+NRkPcLVOGxsfH4+5edhxBZGQk48aNw8vLi2XLljF16lRWrlxZrq9JkyYRGRlJZGQkzs7O+hxSrbDqHAxAzm2X0Xq6a+eQ+P7Y9xQU6/cH2drchAPn0xj51R6+2xVbozjLsGimLfux6zPD9SmEaFD0ujAeHR3N0aNHS9/369ePTp066dRHaGgoMTExxMbG4uHhwZIlS/jtt9/KtImNvfkHb+LEiQwbNoyRI0fqE3K9YOHtjVVwMDl795EXfRpLv3al61rbtebCtQtcyb6Cp51nFb1UbO4jnYm+nMmTP0aQlV9ouKAf+gXmdII8efRZiKZKrzOKzp07s2/fvtL3+/fvp1evXjr1YWpqyty5cxk8eDD+/v6MHTuWgIAA5s+fz/z58/UJq0Gw8PUBIO2H78ssfzboWQD+OP0H5zJ0v9fgZm9FXz8X1BU8jVYjzb3Awh5iNsLSx7Sv8C8Nuw8hRL2m0xlFYGAgKpWKwsJCfvrpJ1q1agXAxYsXyz2tVB333Xcf9913X5llld24/vHHH3Xuvz5y+/BDsvftJ2tPOCXZ2aVPP7W0aYmJyoSFxxdyJecKH/f5uI4jvYX//ZAQCcnR2nIf53dDz+frOiohRC3RKVE09WlODcXEzo7CuDhSFnyLy79eAtDWfBofwai/RpGcq/8N+aIShW93xvLyQD/UagOdXYz86ubva1+G4ysM068QokHQKVHcOv1peno6MTExZZ52un16VFExz28XENOzFyWZZa/7m5mYkV2YTcRl/ceMWJmZkFtYTEJGLp6O1jUNtWI5afB5IJhawIM/QsuOxtmPEKJe0OsexXfffUefPn0YPHgw7777LoMHD+a9994zcGiNl6mjI5iZkf7bYopSy9Z+6uLSBbVK/+EtH47U/tHeey6VbdFJbItO4ljC1RrFW0bwI9qXW5B2boukE4brWwhRL+n11NOcOXOIiIggLCyMbdu2cerUKd59911Dx9ao2Q+/n6vL/6QoNRXTFi1Kl3vYemCmNtO7X1tL7T/pa8uiSpepVHDknUHYW+nfbymPrtpX6lntCG4hRKOn11dXS0tLLC0tAcjPz6d9+/ZER0cbNLDGrtld2tIZl954E6W4+A6tq2+gvytrX+jNiqk9WTG1J8/e7Y2iQH6R4fYBwI2znj+fgffs4QMnOPu3YfchhKgX9Dqj0Gg0ZGRkMHLkSAYOHEjz5s3LDZYTVbMK1F4iyjtxguKMjJtnFQrkF+dzOfsyLW1a6tyvWq0iwN2+9P2JS9r7IH0/3Y5apWJQgCv/Gxtc4/hp7gX3fgo5qVCQBXvnQlostK1510KI+kWvRLFihfapl/fee49+/fpx9epVhgwZYtDAGjszDw9c33mbKx98WGb5jeRw7uo5vRLF7QZ2cCUuLZfC4hL+PpXE0biMGvcJaK9ndZ+k/T3zijZR5GdCVlLZdpYOYGpumH0KIeqEXokiLy+Pr7/+mt27d6NSqejduzclJTWoLyRK+bfwB2D+0fkk5SQx0mdkjfpzsbXkjXvbA3D5Wh57zqTw2cablwld7CyYENa6wrLx1XZjXost72pft3LvApO26d+3EKLO6ZUoHnvsMWxtbXn+ee2gq8WLFzNhwgT++OMPgwbXVOQePUqzvn1RqdW0tmtNsHMw0WnRLD61uMaJ4lZ+rrZsOHaZeTu0kxyVKAqKAkMCWuJiZ6l/x9aO8NCvkHW57PKopZBxsQYRCyHqgzqr9STApFkzAOKnPkerhT9g06MHjpaO/Hzfzzy39TlSclMMur8XBvjywgDf0vdLDlzkjT//obiCSr468x9Wftmlo9qCgkknwbm99nKVEKLBqbNaTwLshg3D/dNPASjJza31/d9ID3vOGGYe73LMrLU3u78Og3/kbFOIhqpOaz01dSq1GnPvNpWuP5Fq3MFsPby1T1rlFhQZZwd3vw7unWHFs5CbYZx9CCGMTmo91VNZBVkApOSm4GTlZJR9NLs+OG/ZoQQOX8woXe7rasuUvgZ4ztXaEXwGan/PuAAFOWBupLIiQgij0enSU+vWrUtfGRkZrF69mtWrV5ORkSF1nmooe+++Mu+HtdVe8197bi35xflG2ae9lRmhXs1Jy84n4kIaERfS2HzyCp9uPGW4nZiag8pE+/jspv8Yrl8hRK3R6x7FnDlzGD9+PElJSSQlJfHoo4/y5ZcyR4E+LNpoLz2V5OaUWe5q7QrAZ5GfsTdxr1H2bWai5o/JPdn1Wv/S18SeXijAgdi00ld2fg0uTVnYwqTt0MxVO85CCNHg6PXU0/fff8/+/fuxuT6Xwuuvv06PHj1KH5cV1ae2tsbE2Ymry5bT4qmnShNHH00ffr73Zyasn0B4YjhBzkE4WjoaPR5rc1MUBcZ+czM5PdK9FR+NCtS/U7cg7Y3tawlwam3l7Zz8wMlH//0IIYxCr0ShKAomJial701MTFAM8YhlE2Ud3JnMzZu5unJV6fwUAG42bpioTFh8ajG25rY839n4ifiJXl4EezpQcv3f8+WlR8nKM8DNbkt7uLBH+6qMkx9MO1DzfQkhDEqvRPHEE0/QvXt3Ro0aBcDKlSt56qmnDBpYU6L58v842d6f7L174ZZE4Wrjyo6HdtB/aX8SshI4mXoSP0e/GpUhvxNLMxN6tL1ZzdbCTM1fRxMZ01XD3e2c9e/4sZVVD77b9hFckZLlQtRHeiWK6dOn07dvX3bv3o2iKCxcuJDOnTsbOrYmJy8qisIrVzBzdS1dZm9hj42ZDWvPrWXtubV80e8LBrQaUGsxPdytFbPXn2J7dFLNEoVVc+2r0vWOUFxQdbJwaAUWzfSPQQihF50ThaIoxMfH06VLF7p06WKMmJok17fe4srMmRQmJGDq7IxKffOsYdG9iziafJS397zNtfxrVfRieJPvbstX285wMTWHTce1JTq8nW3wcbE17I7MLLUlQOb1qLxNm7vh8b8Mu18hxB3pnChUKhUjR47k4MGDxoinyTJtob1RfeGR8TSfMIGWb/27dF0b+zZYm2rHH6w+t5pRvqNqNTY7SzO2nkpi6yltZVgPByv2vNHfsDvp95Y2EVRm9+eQm27YfQohqkWvS09hYWFEREQQGhpq6HiarGb9+uH+8WyufPoZRSnJ5da72mgvR9mb25dbZ2x/TevFpavaudHnbT/L3nNGKPlh4wQBIytff3SJ9qkpIUSt0+uu6LZt2wgLC6Nt27YEBQURGBhIUFCQoWNrUtRWVtiPGIHaxprM9RsqbOPX3I8tF7fQ+afOrD1XxWOmBtaimQUdPezp6GGPo405WXlFTP/9CG/+GUXStbzaCUKlgivH4VPfsq9Fw2tn/0I0YXqdUaxfv97QcYjrTOzsKVRXnL+nh0zn0JVDfBP1Deevna/dwK4L8WrOjtPJ7D2XyqWreXRv04KRnT2Mv+OwKdpBe7dKOAjndxl/30I0cXolCldX13ITF02ZMsXQsTVJNt27kX/6dIXrerr3pKd7T76J+oa/zvzFc8HP1XJ0MCLYgxHBHsSmZNPvs+21t+M2fbSvW/09A64cq70YhGiiZOKiBshUbYqtuYGfOtLRjZkl/rPyGDPWnsDMRM2XD3cmxMv4o8fLUEpg5dTyy03MoM9rYF8LZztCNHIycVE9oygKSn4+hZcvY9ay4jmz7/K4i21x25h/dD6TO02u5Qi1Wjla80J/H1KzC8gtLObPQwmcupxZu4nCIwQcWkPszrLLS4og85K2xHnXibUXjxCNlF6J4sbERWFhYYBMXGRIZi3dAMg/e7bSRDHefzzb4rZxOOlwbYZWhlqtYvogPwCSMvP481ACiRm5HEu4Wqadu4MVjjbmxgnCb4j2dbtrifA/f+PsU4gmSK9EsX///nITF/n7+5dObBQVFWXQIJsSywDtBFCX330P+9GjcJ5a/rJKd7fuBDkFEZ4YjqIoqOp4ilELExNUKvh6+1m+3n62zDpfl2Zsnl7F+AhjSjwC9luhbX+ZhlWIGtArUWzYUPHjm6LmLNr5YT9iBFl79pC19e8KEwWAuYn2W3pBSQEWJha1GWI59tZmrJjai+TMsvNm/LT3PKev1EFpcXMbUJvBwYXa17O7tBVshRB60StRyCRFxmPSzAb3j2cT9+xkilJSKm3Xy6MXkVciGfPXGBwsHPhm4DfYmNnUYqRlBXs6lFu29eQVDl5I57+bopnQozUutpa1E4ylPfzrGESvgzX/gqJaGushRCNlvDKk1bBhwwb8/Pzw8fFh9uzZ5db/+uuvBAUFERQURM+ePcvcQG8Sqijd3s+zH0O9h+Jg4cDR5KNczr5ci4FVT/uWthSVKHz59xk2Hqvl+GxbaosIgrbYoJTBF0JvdZYoiouLee6551i/fj0nTpxg8eLFnDhRtnJomzZt2LFjB1FRUbz99ttMmjSpjqKtfUpREXknTlCSm1vh+rYObZl912we7fAoAAuPLazN8KplYq827HtTW+m2pC7+TquvnzD/OBRWv1AHAQjROBgsUVy+rNs3xgMHDuDj44O3tzfm5uaMGzeOVatWlWnTs2dPmjfXlqYOCwsjPj7eUOHWe2aeGgAKLlyosl1Xl64AxF6LNXpMNfHl32dYGhlXuzv1DIMhH2vPLNLq9+cjRH1msESh68RFCQkJeHp6lr7XaDQkJFRe9O3777/n3nvv1Tu+hsaqUzAAsSNHkRMZWWk7Z2tnerj1ICo5ij9j/qyl6KqvubUZz9zVhpyCIsLPVH7PxSjMLCFsMth73rmtEKJSet3MrsjatboVqato6tTKHvPctm0b33//Pbt3765w/YIFC1iwYAEAycnlK682RHZDBlOcmkLSZ/8la8dOLNq3x6RZxZP2jPQZyd5Le9ket53RvqNrN9A7UKlUvDW0A5tOXCElq4C9ZyuvPOvlZI2bvZXhg1AUSImGogIwNdKYDiEaMb3OKF5//fVqLauKRqMhLu7mpYj4+Hjc3d3LtYuKiuLpp59m1apVtGjRotx6gEmTJhEZGUlkZCTOzjWYha0eUVtZYXf9DCr1229JW/hjpW3v874Pv+Z+HE46zFMbn+Lglfo3V4iNuSm7z6Tw8Lf7Kn1N/CHCODs3s4LsJFj/mnH6F6KR0+uMYvPmzXz88cdllq1fv77csqqEhoYSExNDbGwsHh4eLFmyhN9++61Mm4sXLzJ69Gh+/vln2rVrp0+oDZqZhwdtt2zh3H33kXv4EBnLl2N3772ora3Ltb2/7f1sj9vOgcsH6JTQia6uXWs/4Cr8MDGU2JTsStfP23GWs0lZxtn5/XPgi44y8ZEQetIpUcybN4+vv/6as2fPlpl/IjMzk549e+q2Y1NT5s6dy+DBgykuLubJJ58kICCA+fPnAzB58mQ++OADUlNTmXp90JmpqSmRVVyvb4zMNR6YurUkO3wv2eF7UVlYYj9saLl2jwc8zuMBj9P5p/o5d3lLe0ta2lc+jmLZwXgOX0zniy0VV84F6N6mBT3aVnxWWSUHT3Dy0307IQQAKqWimwWVuHr1Kunp6bz55ptlxj3Y2tri6FjLVUMrERIS0uiSSUl+PvnR0Zwf+xDNH3mYlu+8U2nbzj915omOT/BCl4b1OOi87Wf5eMOpKtsEaez5a1pv/XYwt5u2jEfAbdPIqk2g8wTtuAshmrCq/nbqdEZhb2+Pvb09o0ePxtHREVtbW2bMmMGhQ4d4++236dy5fn6bbejUFhZYtG0LQOGVpCrbFilFfPvPtzwT9AxWpka4MWwkU/q2ZfLd3pWuf3pRJJdrMpueczs4uRq2zyq/ztQKek7Tv28hGjm97lF8+OGHPPjgg+zevZuNGzfyyiuvMHnyZPbv32/o+MR1ahsbLHx9yY2M5MKjE3B+6UWsQ0LKtevQogMnUk+QmpuKxlZTB5Hqr6rihioVXEzN4elFZW94t3O15bUh7e/c+difyy8ryIZZHqAU6xqqEE2KXk89mZiYANpHYqdMmcKIESMoKCgwaGCiPPvRo7Ho4E9OZCTZ4eEVtnnUXztS+/4V99P55850/rkzg5cNJq+B1zvq396VVi2suXQ1r/R1JO4q83ecvfPGoM00Fb2EEHek1xmFh4cHzz77LJs3b+b1118nPz+fkpISQ8cmbtPiiYm0eGIiJ9v7k/HHMpxfKH8f4i6Pu5jaaSoFJdrEfTL1JHsS93Am4wyOlo60sGpR59Vm9fFI91Y80r1VmWX/3RTNV9vO1FFEQjQdeiWKpUuXsmHDBl555RUcHBy4dOkSn376qaFjE5VQ29lRlJxc4Sx4DpYOTAm+OX/5ipgV7Encw8NrHwYgxDWEhUPqX10ofajQ1pAK/mBTuXVmJmrmP9qFrq3rx0MWQjRkeiUKKysrsrOzWbx4Me+88w6FhYU4ODgYODRRGdsBA7i6YgWXP5yB51dzq2w72GswZiZmFBYXsjR6KQlZCURevvlkg6OlI94Old9Ers9GdvYgM7+IktsqDmYXFLPsYDynr2RVL1GEz4Wjv+sfiEoN97wLvgP170OIekyvRDF16lTUajV///0377zzDra2towZM4aICCONrBVltHz/PTK3biVr69Y7trU2s2aY9zAADicd5tiZYzyx8YnS9WqVmp0P7cTewt5o8RqLt3Mz3r0/oNzyy1fzWHYwnotpOZxJysTHxbbiDsysocc0SD9fs0BOrYXzuyRRiEZL76lQDx06VPo4bPPmzeVmdi1Sm5tj5uZG/rVrlGRno7KyQqW+83MJr4a+ylDvm4P1dsbv5KcTP5Gam9ogE0VlLM3UqFTasRnztp9l87/64OtaQbJQqWDwzJrvcIZrzfsQoh7T66knMzMziouLSx9nTE5ORl2NP1TCcOyH3w9AdNcQTncPoyj9zuUpbM1t6e7WvfQV0EL7bXzfpX1GjbW2OVibs3pab14eqC37ci2vqI4jEqJh0+uM4oUXXmDUqFEkJSXx1ltvsWzZMmbMmGHo2EQV7IcPRykuIe/kCTLXbyA7PBz7oeVLe1Slu1t3ANaeW8uJ1JuTRvk4+DCx40RDhlvrOnrYk5qtPcv9YU8suQXF9PZ1quOohGiY9EoU48ePp2vXrmzduhVFUVi5ciX+/v6Gjk1UwdTZGadJz5ATGUnm+g2kLfxR50RhZ25HJ+dOJOUkkZyrLc+eWZDJ6nOrG3yiAGjlaI2rnQXr/7lEWlaBcRPFlRNw5Lc7t6spTTdw8jH+foS4hU61nhqCxljr6U5ixz5EXlQULT94n+Zjx9aor7mH5/JN1Dc4WDiULpsUNIkJHSbUMMq6M/abvahVsGRSD+Ps4POOcLWWZu9rczc8/lft7Es0KQar9STqpxYTHydh+stcXfUXlu3aYdmhAypz/SboGeo9lKzCLIpLtGUtVp9bzeJTi3nI7yHMTRropD8K7ItN40hcBsGeDobvf+peyKl8QiaDWf4MFMtDI6L2SaJoBOzuu4+0X34l9+BBzo97GOeXp+P0zDN69dXGvg1vdHuj9P2R5COcSjvFgqgFTOvcMAvnBWnsOXA+je92nWPuI10MvwMLW+3L2MwsobjQ+PsR4jaSKBoJj8//R37MGeImTaIkq/IJgnQ1p98cBi8fzN9xf5ORnwFAR6eOjPQZabB9GNt/hnVg++lkGtdFViFqjySKRsLM1RUzV1coKSH3yBGD9evezJ27NXcTlRzFpvObyC7MZvOFzQ0qUdyw52wKo7/eU2aZqVrNW0P96WSMS1LGcDUBdn9+53Y2zhA8XgofCoOQRNEI5ezfT1FyMqYGmj987oCbZUJm7JvB0uil3Lv8XixMLPjk7k9o17z+T1M7LtSTHaeTyywrKlbYey6VfedSG0aicPSG2J2w5b3qtW9zt3Z2PyFqSBJFI+Py6iskffoZZ4fdXzpa29TFBa8/lqLW8wb3re5vez85hTlkFmayPW47O+N3oigKZiZmtLFrU+WcEnXp6bu8efqusjWtcgqK6PDOxjqKSA/DvoAhs+/YjKilsPoFKJGBhsIwJFE0MvYjRlCUlIxSqH06Jv90DDmRkVx+511UVpbYdO+O3ZAhevffybkTnZw7cTn7MtvjtjPn0BzmHJoDwKd9PmVIG/37FnegUoFZNWYtbKhPp4l6SxJFI2Pq5ITrmzefWso5eJCE6S+TtXMnJZmZZO3YgYn9zbpOJo4tsPTT/dJRS5uWLBy8kIz8DK4VXOPd8HdZenoph5MOl2kX5BxUpr5UfbTicAJH4zP03l6lUvFsH2+CNA4Gi0mI+kQSRSNn3bUrvju2A5D4n/9wddlyLj7x5M0GKhXt9oZjokeZ+JCW2qlY84vz+fnEz0SnRROdFl26Prcoly0XttDPsx/WZtY1OQyjsDQ1YUB7Fy6m5RBzJUvvfmKSstA0t6p/iSLxEGRerusohLGZWoBbMBix3p4kiibE9Y03cRg5svR95t/bSPvhB7J278F+mP7f+i1MLFgxYkW55R/u/ZClp5cS9lsYc/rNoV+rfnrvwxjUahXfTwytcT9+/1lvgGgMyNxG+3PZk1W3E43HI39Au0FG614SRRNi0swG65CQm++bNyfthx9I+t9/Sfvxx9LlFr6+uM/6qMb7e6LjE7hYuzD3yFwWHl9Ic8vmBLsE17hfcQd+98ET66Eov64jEcZ2NR7+mgb514y6G0kUTZh5mzY4PPgghUlXSpcVxJ7n6ooVlGTfHLTn8OCDNLurt879a2w1PBn4JNvitnE0+ShfHfmKZwIrHzHeyq4VLW1aVrq+vlKrVHy78xwL95w3+r4eDvXk/REdq25kYgqtexo9FlEPJJ+uld1IomjCVGo1bh9+UGZZ1u49JH36KQWx5wDIjz1P5qZNtFm1Sq+b3mZqM5YMW8KDqx9k36V9Vc594ePgU+ElrPruo9Edib6s/z2O6lp9NJETl4z7zVGIikiiEGU0692LZr17lb5PfPPfXF2xgtgRI3CcOBEAy4AA7O8fplO/Xw34igvXLlS6/vtj37M3cS99f+9baZshbYaUqUNVX4zqrKmV/UTFZ1BYXFIr+xLiVpIoRJXcZ31E0ZUr5EREkLF0KSU5OQAUnD+PysKC5g+Pw8T2zgXxXKxdcLF2qXS9uYk57jbula7fnbCbNefWkFeUV2a5t703jwU8Vs2jafiOxl/lrk/+vmM7RxsLljwThpW5SS1EJepcbjpcS9SWbjExM3j3kijEHbX64fvS39N++ZUrM2aQ8tVXAFxdsQK7e+9FZW6Gw0MPYdq8uV77uDGQrzI/HPuBX0/8yq74XaXLsgqzyC3KxdLUslx7M7UZg70G18vHcvX1RK82tLQrf6y3u5CWw8EL6SRn5tOqReM5flGBG0lh3SvaV/thMO5Xg+9GJi4Seim+do2zQ+6lOC2tdJna2hr7kSPKtLPq1An7ESNu39wgfj/1OzP2Vz4F74QOE+jn2Y8Q15B6W1rEGJYfjOflP46y89V+kigaO0WBk6shNw32zQPzZvDMVr26qupvpyQKUWMlBQWcf3AshZcuoTK5eamjOD0dAJtevXB5/TUs2xm+eGBKbgq3/yecXZjNyFUjKVa0ky9NC57GKN9RVV76akxuJApPRyvMTIw3CEvUDxamJnz5cDA+Gx+HvKtGSRRy6UnUmNrcHO9VK8stz9y2jeTPvyB7zx5ih4/Apnfv0hviFTH38sJc46HTvp2sys+D7Ywz60avY2f8Tmbun8ncI3P5+eTPrB+9HlvzWphgqI71aNuC0V08KCiSG9+NXVZ+EdujkzlxKRNjzqQuiUIYjW2/ftj260faz79wZeZMsnfvJnv37iq3cXjwgUrX2fTqjd2QwdXat3szdx7ye4hubt14dvOzXM6+zD1/3FPuqSlTtSl9Pfs2qgTi7mDF/8YG13UYohacScpie/QOo++nTi89bdiwgRdffJHi4mKefvpp3nij7P/EiqLw4osvsm7dOqytrfnxxx/p0qXqqSzl0lP9VJSSQsHFuErXZ/zxB1k7d5a5dFVm+6Qk7S+33mu4/p+uy2uvlWuvMjPDfuQITGxtycjL4N4/7yWrsOKxDuP9xzOg1YBqHkl5KlQEOgdiYWKhdx9C6ONMUhb3/G8H3s42/C//A9wtC3CZXvWXscrUy3sUxcXFtGvXjs2bN6PRaAgNDWXx4sV06NChtM26dev48ssvWbduHfv37+fFF19k//79VfYriaJxyv3nGFnbyj4WenXNWgovXrzjts0GDMAitAuFw/ujumVOjqKSIkatGkVBSUGN43O3cefJjsavrdTZtXODmChK1I68wmJeXnqU9JwCnkt4DSfTPPz+c0CvvurlPYoDBw7g4+ODt7d2Mplx48axatWqMoli1apVPPbYY6hUKsLCwsjIyODSpUu4ubnVVdiijlgFdsQqsGzpCucXXqAkJwel5PbvOgoZvy8l99g/ZG7YSNbWrWRt3QqzP8XM4+Y9EEUp4Zd8C4o8W1PS0Vfv2NadWU2qbTxrj35w58Y1tNxcRTOXSu7jXD/ZUlC4mn+Nse0exNzUvOzKCiiVrdL3SbHKtquyu/odn6LPvqraRo/PorLP4S5f6OIZxoUfTcgrMs7TfXWWKBISEvD0vDlNo0ajKXe2UFGbhIQESRSilNq64sc/Wzx1/dv951CYmEjKvHkoBYVl2ijFxVxbswbT1HT4R/+aOSOLFaA2T8zvfBal9Z1RoxD1RzGgAfb6N6fy0Uj6q7NEUdEVr9ufda9OG4AFCxawYMECAJKTk8utF02bmbs7bh9+WOE6j88+rXH/hVeSKLpi/HkfCq9cQckvqPgL523/q+QW5VCsXH/qqaqry/pcedanP31jqHSVfv1VeqW9qo+hys9In+PVY1+V9qcQf/owhfZWZO2PxKnvXVV0rr86SxQajYa4uJs3N+Pj43F3d9e5DcCkSZOYNGkSoL3OJkRtMnN1wczV+GM0qjEJain7OzcRjURtVBqrs9E4oaGhxMTEEBsbS0FBAUuWLGH48OFl2gwfPpyffvoJRVHYt28f9vb2ctlJCCFqWZ2dUZiamjJ37lwGDx5McXExTz75JAEBAcyfPx+AyZMnc99997Fu3Tp8fHywtrZm4cKFdRWuEEI0WVLCQwghRJV/O6UQjBBCiCpJohBCCFElSRRCCCGqJIlCCCFElSRRCCGEqFKje+rJyckJLy8vvbZNTk7G2dnZsAHVc3LMTYMcc9NQk2M+f/48KSkpFa5rdImiJprio7VyzE2DHHPTYKxjlktPQgghqiSJQgghRJUkUdziRmHBpkSOuWmQY24ajHXMco9CCCFEleSMQgghRJUkUVy3YcMG/Pz88PHxYfbs2XUdjsHExcXRr18//P39CQgIYM6cOQCkpaUxcOBAfH19GThwIOnp6aXbzJo1Cx8fH/z8/Ni4cWNdhV4jxcXFdO7cmWHDhgGN/3gBMjIyeOCBB2jfvj3+/v7s3bu3UR/3559/TkBAAB07duThhx8mLy+vUR7vk08+iYuLCx073pwKWJ/jPHjwIIGBgfj4+PDCCy9UPolTRRShFBUVKd7e3srZs2eV/Px8JSgoSDl+/Hhdh2UQiYmJysGDBxVFUZRr164pvr6+yvHjx5VXX31VmTVrlqIoijJr1izltddeUxRFUY4fP64EBQUpeXl5yrlz5xRvb2+lqKiozuLX13//+1/l4YcfVoYOHaooitLoj1dRFOWxxx5Tvv32W0VRFCU/P19JT09vtMcdHx+veHl5KTk5OYqiKMqDDz6oLFy4sFEe744dO5SDBw8qAQEBpcv0Oc7Q0FAlPDxcKSkpUYYMGaKsW7eu2jFIolAUJTw8XBk0aFDp+48++kj56KOP6jAi4xk+fLiyadMmpV27dkpiYqKiKNpk0q5dO0VRyh/7oEGDlPDw8DqJVV9xcXFK//79la1bt5YmisZ8vIqiKFevXlW8vLyUkpKSMssb63HHx8crGo1GSU1NVQoLC5WhQ4cqGzdubLTHGxsbWyZR6HqciYmJip+fX+ny3377TZk0aVK19y+XnoCEhAQ8PT1L32s0GhISEuowIuM4f/48hw8fpnv37ly5cqV0tkA3NzeSkpKAxvFZvPTSS3zyySeo1Tf/827Mxwtw7tw5nJ2deeKJJ+jcuTNPP/002dnZjfa4PTw8eOWVV2jVqhVubm7Y29szaNCgRnu8t9P1OBMSEtBoNOWWV5ckCiqecF2lqmgG+4YrKyuLMWPG8MUXX2BnZ1dpu4b+WaxZswYXFxe6du1arfYN/XhvKCoq4tChQ0yZMoXDhw9jY2NT5b22hn7c6enprFq1itjYWBITE8nOzuaXX36ptH1DP97qquw4a3r8kijQZte4uLjS9/Hx8bi7u9dhRIZVWFjImDFjGD9+PKNHjwbA1dWVS5cuAXDp0iVcXFyAhv9Z7Nmzh7/++gsvLy/GjRvH33//zaOPPtpoj/cGjUaDRqOhe/fuADzwwAMcOnSo0R73li1baNOmDc7OzpiZmTF69GjCw8Mb7fHeTtfj1Gg0xMfHl1teXZIogNDQUGJiYoiNjaWgoIAlS5YwfPjwug7LIBRF4amnnsLf35/p06eXLh8+fDiLFi0CYNGiRYwYMaJ0+ZIlS8jPzyc2NpaYmBi6detWJ7HrY9asWcTHx3P+/HmWLFlC//79+eWXXxrt8d7QsmVLPD09iY6OBmDr1q106NCh0R53q1at2LdvHzk5OSiKwtatW/H392+0x3s7XY/Tzc0NW1tb9u3bh6Io/PTTT6XbVIt+t1Yan7Vr1yq+vr6Kt7e3MmPGjLoOx2B27dqlAEpgYKDSqVMnpVOnTsratWuVlJQUpX///oqPj4/Sv39/JTU1tXSbGTNmKN7e3kq7du10ejKivtm2bVvpzeymcLyHDx9WunbtqgQGBiojRoxQ0tLSGvVxv/POO4qfn58SEBCgPProo0peXl6jPN5x48YpLVu2VExNTRUPDw/lu+++0+s4IyIilICAAMXb21t57rnnyj34UBUZmS2EEKJKculJCCFElSRRCCGEqJIkCiGEEFWSRCGEEKJKkiiEEEJUSRKFEEKIKkmiEEIIUSVJFELoKCMjg6+//rr0fc+ePY2yn/j4eH7//Xej9C2ELiRRCKGj2xNFeHi4UfazdetWDh06ZJS+hdCFJAohdPTGG29w9uxZgoODefXVV2nWrBmgLePevn17nn76aTp27Mj48ePZsmULvXr1wtfXlwMHDpT28csvv9CtWzeCg4N59tlnKS4uLrOP3bt3M336dJYtW0ZwcDCxsbG1eoxClGHAkiRCNAm3TyJjY2NTutzExESJiopSiouLlS5duihPPPGEUlJSoqxcuVIZMWKEoiiKcuLECWXYsGFKQUGBoiiKMmXKFGXRokXl9jN48GDln3/+Mf4BCXEHpnWdqIRoTNq0aUNgYCAAAQEBDBgwAJVKRWBgIOfPnwe0l5QOHjxIaGgoALm5uaVlom8VHR2Nn59frcUuRGUkUQhhQBYWFqW/q9Xq0vdqtZqioiJAW/r98ccfZ9asWZX2k5qair29PWZmZsYNWIhqkHsUQujI1taWzMxMvbcfMGAAy5YtK52+Mi0tjQsXLpRpExsb26An1hGNiyQKIXTUokULevXqRceOHXn11Vd13r5Dhw7MmDGDQYMGERQUxMCBA0tnK7uhffv2pKSk0LFjR6M9VSVEdcl8FEIIIaokZxRCCCGqJIlCCCFElSRRCCGEqJIkCiGEEFWSRCGEEKJKkiiEEEJUSRKFEEKIKkmiEEIIUaX/B0vUoKRS+X+tAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "pred_surv = estimator.predict_survival_function(x_new)\n",
    "time_points = np.arange(1, 1000)\n",
    "for i, surv_func in enumerate(pred_surv):\n",
    "    plt.step(time_points, surv_func(time_points), where=\"post\", label=f\"Sample {i + 1}\")\n",
    "plt.ylabel(r\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measuring the Performance of Survival Models\n",
    "\n",
    "Once we fit a survival model, we usually want to assess how well a model can actually predict survival. Our test data is usually subject to censoring too, therefore metrics like root mean squared error or correlation are unsuitable. Instead, we use generalization of the area under the receiver operating characteristic (ROC) curve called [Harrell's concordance index](https://pdfs.semanticscholar.org/7705/392f1068c76669de750c6d0da8144da3304d.pdf) or c-index.\n",
    "\n",
    "The interpretation is identical to the traditional area under the [ROC curve](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) metric for binary classification:\n",
    "- a value of 0.5 denotes a random model,\n",
    "- a value of 1.0 denotes a perfect model,\n",
    "- a value of 0.0 denotes a perfectly wrong model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.73626'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.metrics import concordance_index_censored\n",
    "\n",
    "prediction = estimator.predict(data_x_numeric)\n",
    "result = concordance_index_censored(data_y[\"Status\"], data_y[\"Survival_in_days\"], prediction)\n",
    "f\"{result[0]:.5f}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or alternatively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.73626'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_index = estimator.score(data_x_numeric, data_y)\n",
    "f\"{c_index:.5f}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model's c-index indicates that the model clearly performs better than random, but is also far from perfect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Selection: Which Variable is Most Predictive?\n",
    "\n",
    "The model above considered all available variables for prediction. Next, we want to investigate which single variable is the best risk predictor. Therefore, we fit a Cox model to each variable individually and record the c-index on the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Karnofsky_score          0.709280\n",
       "Celltype=smallcell       0.572581\n",
       "Celltype=large           0.561620\n",
       "Celltype=squamous        0.550545\n",
       "Treatment=test           0.525386\n",
       "Age_in_years             0.515107\n",
       "Months_from_Diagnosis    0.509030\n",
       "Prior_therapy=yes        0.494434\n",
       "dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def fit_and_score_features(X, y):\n",
    "    n_features = X.shape[1]\n",
    "    scores = np.empty(n_features)\n",
    "    m = CoxPHSurvivalAnalysis()\n",
    "    for j in range(n_features):\n",
    "        Xj = X[:, j : j + 1]\n",
    "        m.fit(Xj, y)\n",
    "        scores[j] = m.score(Xj, y)\n",
    "    return scores\n",
    "\n",
    "\n",
    "scores = fit_and_score_features(data_x_numeric.values, data_y)\n",
    "pd.Series(scores, index=data_x_numeric.columns).sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Karnofsky_score` is the best variable, whereas `Months_from_Diagnosis` and `Prior_therapy='yes'` have almost no predictive power on their own.\n",
    "\n",
    "Next, we want to build a parsimonious model by excluding irrelevant features. We could use the ranking from above, but would need to determine what the optimal cut-off should be. Luckily, scikit-learn has built-in support for performing grid search.\n",
    "\n",
    "First, we create a pipeline that puts all the parts together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "pipe = Pipeline(\n",
    "    [\n",
    "        (\"encode\", OneHotEncoder()),\n",
    "        (\"select\", SelectKBest(fit_and_score_features, k=3)),\n",
    "        (\"model\", CoxPHSurvivalAnalysis()),\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to define the range of parameters we want to explore during grid search. Here, we want to optimize the parameter `k` of the `SelectKBest` class and allow `k` to vary from 1 feature to all 8 features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>param_select__k</th>\n",
       "      <th>params</th>\n",
       "      <th>split0_test_score</th>\n",
       "      <th>split1_test_score</th>\n",
       "      <th>split2_test_score</th>\n",
       "      <th>mean_test_score</th>\n",
       "      <th>std_test_score</th>\n",
       "      <th>rank_test_score</th>\n",
       "      <th>split0_train_score</th>\n",
       "      <th>split1_train_score</th>\n",
       "      <th>split2_train_score</th>\n",
       "      <th>mean_train_score</th>\n",
       "      <th>std_train_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>{'select__k': 5}</td>\n",
       "      <td>0.716093</td>\n",
       "      <td>0.719862</td>\n",
       "      <td>0.716685</td>\n",
       "      <td>0.717547</td>\n",
       "      <td>0.001655</td>\n",
       "      <td>1</td>\n",
       "      <td>0.732087</td>\n",
       "      <td>0.742432</td>\n",
       "      <td>0.731710</td>\n",
       "      <td>0.735410</td>\n",
       "      <td>0.004968</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>{'select__k': 4}</td>\n",
       "      <td>0.697368</td>\n",
       "      <td>0.722332</td>\n",
       "      <td>0.727324</td>\n",
       "      <td>0.715675</td>\n",
       "      <td>0.013104</td>\n",
       "      <td>2</td>\n",
       "      <td>0.732477</td>\n",
       "      <td>0.743090</td>\n",
       "      <td>0.727138</td>\n",
       "      <td>0.734235</td>\n",
       "      <td>0.006630</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>8</td>\n",
       "      <td>{'select__k': 8}</td>\n",
       "      <td>0.706478</td>\n",
       "      <td>0.723320</td>\n",
       "      <td>0.716685</td>\n",
       "      <td>0.715494</td>\n",
       "      <td>0.006927</td>\n",
       "      <td>3</td>\n",
       "      <td>0.739356</td>\n",
       "      <td>0.746249</td>\n",
       "      <td>0.737519</td>\n",
       "      <td>0.741041</td>\n",
       "      <td>0.003758</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6</td>\n",
       "      <td>{'select__k': 6}</td>\n",
       "      <td>0.704453</td>\n",
       "      <td>0.719368</td>\n",
       "      <td>0.716685</td>\n",
       "      <td>0.713502</td>\n",
       "      <td>0.006491</td>\n",
       "      <td>4</td>\n",
       "      <td>0.735722</td>\n",
       "      <td>0.747565</td>\n",
       "      <td>0.731710</td>\n",
       "      <td>0.738332</td>\n",
       "      <td>0.006731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7</td>\n",
       "      <td>{'select__k': 7}</td>\n",
       "      <td>0.700405</td>\n",
       "      <td>0.719368</td>\n",
       "      <td>0.720045</td>\n",
       "      <td>0.713272</td>\n",
       "      <td>0.009103</td>\n",
       "      <td>5</td>\n",
       "      <td>0.741173</td>\n",
       "      <td>0.742564</td>\n",
       "      <td>0.728621</td>\n",
       "      <td>0.737453</td>\n",
       "      <td>0.006271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>{'select__k': 2}</td>\n",
       "      <td>0.699393</td>\n",
       "      <td>0.717885</td>\n",
       "      <td>0.718365</td>\n",
       "      <td>0.711881</td>\n",
       "      <td>0.008833</td>\n",
       "      <td>6</td>\n",
       "      <td>0.732087</td>\n",
       "      <td>0.727428</td>\n",
       "      <td>0.714409</td>\n",
       "      <td>0.724642</td>\n",
       "      <td>0.007481</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>{'select__k': 1}</td>\n",
       "      <td>0.698887</td>\n",
       "      <td>0.707510</td>\n",
       "      <td>0.712206</td>\n",
       "      <td>0.706201</td>\n",
       "      <td>0.005516</td>\n",
       "      <td>7</td>\n",
       "      <td>0.710670</td>\n",
       "      <td>0.714793</td>\n",
       "      <td>0.700445</td>\n",
       "      <td>0.708636</td>\n",
       "      <td>0.006032</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>{'select__k': 3}</td>\n",
       "      <td>0.708502</td>\n",
       "      <td>0.714427</td>\n",
       "      <td>0.694849</td>\n",
       "      <td>0.705926</td>\n",
       "      <td>0.008198</td>\n",
       "      <td>8</td>\n",
       "      <td>0.734034</td>\n",
       "      <td>0.722559</td>\n",
       "      <td>0.716634</td>\n",
       "      <td>0.724409</td>\n",
       "      <td>0.007223</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   param_select__k            params  split0_test_score  split1_test_score  \\\n",
       "4                5  {'select__k': 5}           0.716093           0.719862   \n",
       "3                4  {'select__k': 4}           0.697368           0.722332   \n",
       "7                8  {'select__k': 8}           0.706478           0.723320   \n",
       "5                6  {'select__k': 6}           0.704453           0.719368   \n",
       "6                7  {'select__k': 7}           0.700405           0.719368   \n",
       "1                2  {'select__k': 2}           0.699393           0.717885   \n",
       "0                1  {'select__k': 1}           0.698887           0.707510   \n",
       "2                3  {'select__k': 3}           0.708502           0.714427   \n",
       "\n",
       "   split2_test_score  mean_test_score  std_test_score  rank_test_score  \\\n",
       "4           0.716685         0.717547        0.001655                1   \n",
       "3           0.727324         0.715675        0.013104                2   \n",
       "7           0.716685         0.715494        0.006927                3   \n",
       "5           0.716685         0.713502        0.006491                4   \n",
       "6           0.720045         0.713272        0.009103                5   \n",
       "1           0.718365         0.711881        0.008833                6   \n",
       "0           0.712206         0.706201        0.005516                7   \n",
       "2           0.694849         0.705926        0.008198                8   \n",
       "\n",
       "   split0_train_score  split1_train_score  split2_train_score  \\\n",
       "4            0.732087            0.742432            0.731710   \n",
       "3            0.732477            0.743090            0.727138   \n",
       "7            0.739356            0.746249            0.737519   \n",
       "5            0.735722            0.747565            0.731710   \n",
       "6            0.741173            0.742564            0.728621   \n",
       "1            0.732087            0.727428            0.714409   \n",
       "0            0.710670            0.714793            0.700445   \n",
       "2            0.734034            0.722559            0.716634   \n",
       "\n",
       "   mean_train_score  std_train_score  \n",
       "4          0.735410         0.004968  \n",
       "3          0.734235         0.006630  \n",
       "7          0.741041         0.003758  \n",
       "5          0.738332         0.006731  \n",
       "6          0.737453         0.006271  \n",
       "1          0.724642         0.007481  \n",
       "0          0.708636         0.006032  \n",
       "2          0.724409         0.007223  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import GridSearchCV, KFold\n",
    "\n",
    "param_grid = {\"select__k\": np.arange(1, data_x_numeric.shape[1] + 1)}\n",
    "cv = KFold(n_splits=3, random_state=1, shuffle=True)\n",
    "gcv = GridSearchCV(pipe, param_grid, return_train_score=True, cv=cv)\n",
    "gcv.fit(data_x, data_y)\n",
    "\n",
    "results = pd.DataFrame(gcv.cv_results_).sort_values(by=\"mean_test_score\", ascending=False)\n",
    "results.loc[:, ~results.columns.str.endswith(\"_time\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results show that it is sufficient to select the 5 most predictive features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Celltype=large       -0.754714\n",
       "Celltype=smallcell   -0.328059\n",
       "Celltype=squamous    -1.147673\n",
       "Karnofsky_score      -0.031112\n",
       "Treatment=test        0.257313\n",
       "dtype: float64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pipe.set_params(**gcv.best_params_)\n",
    "pipe.fit(data_x, data_y)\n",
    "\n",
    "encoder, transformer, final_estimator = (s[1] for s in pipe.steps)\n",
    "pd.Series(final_estimator.coef_, index=encoder.encoded_columns_[transformer.get_support()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's next?\n",
    "\n",
    "Cox's proportional hazards model is by far the most popular survival model, because once trained, it is easy to interpret. However, if prediction performance is the main objective, more sophisticated, non-linear or ensemble models might lead to better results.\n",
    "Before you dive deeper into the various survival models, it is highly recommended reading [this notebook](evaluating-survival-models.ipynb) for getting a better understanding on how to evaluate survival models.\n",
    "The [User Guide](https://scikit-survival.readthedocs.io/en/latest/user_guide/index.html) is a good starting point to learn more about various models implemented in **scikit-survival**, and the [API reference](https://scikit-survival.readthedocs.io/en/latest/api/index.html) contains a full list of available classes and functions and their parameters. In addition, you can use any unsupervised pre-processing method available with scikit-learn, for instance, you could perform dimensionality reduction using [Non-Negative Matrix Factorization (NMF)](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.NMF.html#sklearn.decomposition.NMF), before training a Cox model."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
